Fitting the ST model for black carbon in Denver, CO
Sheena Martenies
01/05/2021
1 Introduction
In this script, we are building the spatiotemporal prediction model for black carbon in the Denver metro area. This model uses data from 5 sampling campaigns (as detailed in 15_Exploring_BC_Data.R). Two models are developed here: one which uses one temporal basis function (Model A) and one that uses two temporal basis functions (Model B).
Updates to this model are in response to reviewer comments on an earlier manuscript version submitted to ES&T.
A summary table at the end of the document summarizes the performance for each of these models. Overall, Model B performs well and allows us to predict BC concentrations at a weekly temporal resolution from 2009 through 2020.
## Loading required package: Matrix2 Developing data sets
2.1 Denver data set
Here I’m reading in the filter data sets and getting them into the right shape for the SpatioTemporal package. Here the ST covariates cover the entire study period (2009 - 2020).
data_name <- "Combined_Filter_Data_AEA.csv"
all_data <- read_csv(here::here("Data", data_name))## Parsed with column specification:
## cols(
## .default = col_double(),
## filter_id = col_character(),
## campaign = col_character(),
## StartDateTimeLocal = col_date(format = ""),
## EndDateTimeLocal = col_date(format = ""),
## sample_week = col_date(format = ""),
## As_ug_m3 = col_logical(),
## Cd_ug_m3 = col_logical(),
## In_ug_m3 = col_logical(),
## Ni_ug_m3 = col_logical(),
## Sn_ug_m3 = col_logical(),
## Te_ug_m3 = col_logical(),
## st_week = col_date(format = ""),
## units_pm = col_character(),
## units_bc = col_character(),
## units_no2 = col_character(),
## units_temp = col_character(),
## units_smoke = col_character(),
## nn3_bc = col_logical(),
## area_bc = col_logical(),
## idw_bc = col_logical()
## )## See spec(...) for full column specifications.#' Select a "calibrated" version of the data
#' For now, go with Deming regression-- accounts for variability in the
#' monitor and the UPAS data and the temporal mismatch in TWAs
all_data$bc_ug_m3 <- all_data$bc_ug_m3_dem
all_data$pm_ug_m3 <- all_data$pm_ug_m3_dem
#' List of sites that failed preliminary screening (see 15_Exploring_BC_Data.R)
#' These are all the sites (by campaign) where BC measurements had a CV > 30%
#' Note that all of these are in Campaign 4
cv_drop <- read_csv(here::here("Data/Dropped_Sites_by_Campaign.csv"))## Parsed with column specification:
## cols(
## site_id = col_double(),
## campaign = col_character(),
## n_filters = col_double(),
## mean = col_double(),
## sd = col_double(),
## cv = col_double(),
## median = col_double(),
## IQR = col_double(),
## range = col_double(),
## cv_gt_50 = col_double(),
## cv_gt_30 = col_double()
## )filter_data <- filter(all_data, is.na(filter_id) | filter_id != "080310027") %>%
filter(is.na(indoor) | indoor == 0) %>%
filter(is.na(is_blank) | is_blank == 0) %>%
#' QA filters
filter(is.na(bc_below_lod) | bc_below_lod == 0) %>%
filter(is.na(negative_pm_mass) | negative_pm_mass == 0) %>%
filter(is.na(potential_contamination) | potential_contamination == 0)
dist_data <- bind_rows(filter(filter_data, is.na(campaign) | campaign != "Campaign4"),
filter(filter_data, campaign == "Campaign4" & !(site_id %in% cv_drop$site_id))) %>%
arrange(st_week)
head(dist_data$sample_week)## [1] NA NA NA NA NA NAtail(dist_data$sample_week)## [1] NA NA NA NA NA NAhead(dist_data$st_week)## [1] "2008-12-29" "2008-12-29" "2008-12-29" "2008-12-29" "2008-12-29"
## [6] "2008-12-29"tail(dist_data$st_week)## [1] "2020-12-28" "2020-12-28" "2020-12-28" "2020-12-28" "2020-12-28"
## [6] "2020-12-28"length(unique(dist_data$filter_id))## [1] 848summary(dist_data$bc_ug_m3)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.51 0.99 1.17 1.27 1.46 2.27 41600quantile(dist_data$bc_ug_m3, probs = c(0.05, 0.95), na.rm = T)## 5% 95%
## 0.8689785 2.0648167sd(dist_data$bc_ug_m3, na.rm = T)## [1] 0.3498975sd(dist_data$bc_ug_m3, na.rm = T) / mean(dist_data$bc_ug_m3, na.rm = T)## [1] 0.2764677#' All of the dropped data should be in Campaign 4
dropped_data <- filter(filter_data, campaign == "Campaign4" & site_id %in% cv_drop$site_id)
length(unique(dropped_data$filter_id))## [1] 97summary(dropped_data$bc_ug_m3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.7752 0.6319 1.4887 1.1935 1.9480 2.1754quantile(dropped_data$bc_ug_m3, probs = c(0.05, 0.95))## 5% 95%
## -0.3137664 2.1225874sd(dropped_data$bc_ug_m3)## [1] 0.9202973sd(dropped_data$bc_ug_m3) / mean(dropped_data$bc_ug_m3)## [1] 0.7711097#' central site data
central_data <- filter(all_data, filter_id == "080310027") %>%
arrange(st_week)
head(central_data$sample_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(central_data$sample_week)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"head(central_data$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(central_data$st_week)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"#' Add prefix to site_ids
unique(dist_data$site_id)## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68dist_data$site_id <- paste0("d_", dist_data$site_id)
unique(dist_data$site_id)## [1] "d_1" "d_2" "d_3" "d_4" "d_5" "d_6" "d_7" "d_8" "d_9" "d_10"
## [11] "d_11" "d_12" "d_13" "d_14" "d_15" "d_16" "d_17" "d_18" "d_19" "d_20"
## [21] "d_21" "d_22" "d_23" "d_24" "d_25" "d_26" "d_27" "d_28" "d_29" "d_30"
## [31] "d_31" "d_32" "d_33" "d_34" "d_35" "d_36" "d_37" "d_38" "d_39" "d_40"
## [41] "d_41" "d_42" "d_43" "d_44" "d_45" "d_46" "d_47" "d_48" "d_49" "d_50"
## [51] "d_51" "d_52" "d_53" "d_54" "d_55" "d_56" "d_57" "d_58" "d_59" "d_60"
## [61] "d_61" "d_62" "d_63" "d_64" "d_65" "d_66" "d_67" "d_68"unique(central_data$site_id)## [1] 16central_data$site_id <- "central"
unique(central_data$site_id)## [1] "central"#' put these data sets back together
lur_data <- bind_rows(dist_data, central_data) %>%
arrange(site_id, st_week)
head(lur_data$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(lur_data$st_week)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"#' Log-transformed BC observations with NA values for the missing dates
#' Note that campaign 5 has duplicate measures, so we need to average them
bc_obs <- lur_data %>%
dplyr::select(site_id, st_week, bc_ug_m3) %>%
group_by(site_id, st_week) %>%
summarize(bc_ug_m3 = mean(bc_ug_m3, na.rm = T)) %>%
mutate(bc_ug_m3 = ifelse(is.nan(bc_ug_m3), NA, bc_ug_m3)) %>%
mutate(log_bc = log(bc_ug_m3)) %>%
dplyr::select(-bc_ug_m3)
#' Pivot to a wide data frame
bc_obs2 <- bc_obs %>%
pivot_wider(id_cols = st_week, names_from = site_id, values_from = log_bc) %>%
# names_from = site_id, values_from = bc_ug_m3) %>%
as.data.frame() %>%
arrange(st_week)
rownames(bc_obs2) <- bc_obs2$st_week
bc_obs2$st_week <- NULL
bc_obs2 <- as.matrix(bc_obs2)
bc_weeks <- rownames(bc_obs2)
tail(bc_weeks)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"#' Spatial covariates (scaled)
load(here::here("Results", "BC_LASSO_Model_5C.Rdata"))
lasso_coef_df <- data.frame(name = log_bc_lasso_coef4@Dimnames[[1]][log_bc_lasso_coef4@i + 1],
coefficient = log_bc_lasso_coef4@x)
covars <- as.character(lasso_coef_df$name[-1])
covars ## [1] "tree_cover_100" "impervious_2500" "low_int_100" "med_int_50"
## [5] "high_int_100" "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells"
## [9] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"covar_fun <- paste("~", paste(covars, collapse = " + "))
covar_fun## [1] "~ tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + high_int_100 + pop_ct_1000 + dist_m_cafos + dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000"bc_sp_cov <- dplyr::select(lur_data, site_id, lon, lat, all_of(covars)) %>%
distinct() %>%
mutate_at(.vars = vars(all_of(covars)),
scale)
bc_sp_cov <- bc_sp_cov %>%
rename(ID = site_id) %>%
mutate(lon2 = lon, lat2 = lat) %>%
st_as_sf(coords = c("lon2", "lat2"), crs = ll_wgs84) %>%
st_transform(crs = albers)
sp_coords <- do.call(rbind, st_geometry(bc_sp_cov)) %>% as_tibble()
names(sp_coords) <- c("x", "y")
bc_sp_cov <- bind_cols(bc_sp_cov, sp_coords) %>%
st_set_geometry(NULL) %>%
as.data.frame()
bc_sp_cov$type <- ifelse(bc_sp_cov$ID == "central", "central", "dist")
bc_sp_cov$type <- as.factor(bc_sp_cov$type)
cor(bc_sp_cov[,covars])## tree_cover_100 impervious_2500 low_int_100 med_int_50
## tree_cover_100 1.000000000 0.05592620 0.227568688 -0.41656599
## impervious_2500 0.055926204 1.00000000 0.046627086 -0.13878441
## low_int_100 0.227568688 0.04662709 1.000000000 -0.11614316
## med_int_50 -0.416565993 -0.13878441 -0.116143165 1.00000000
## high_int_100 -0.335160074 0.48006634 -0.444088483 0.11084350
## pop_ct_1000 0.232959389 0.43482586 0.457037018 -0.02240302
## dist_m_cafos 0.262782544 -0.02217726 -0.161007912 -0.09924381
## dist_m_og_wells 0.304301808 -0.11331048 0.008034836 -0.21284246
## dist_m_wwtp 0.140092353 -0.49234185 -0.075238446 0.02725541
## aadt_100 -0.176079221 0.40911376 -0.255918459 0.06413761
## aadt_250 -0.007357196 0.42742607 -0.261080670 0.03197137
## aadt_1000 0.085749987 0.43050282 -0.112576947 -0.19807436
## high_int_100 pop_ct_1000 dist_m_cafos dist_m_og_wells
## tree_cover_100 -0.33516007 0.23295939 0.26278254 0.304301808
## impervious_2500 0.48006634 0.43482586 -0.02217726 -0.113310476
## low_int_100 -0.44408848 0.45703702 -0.16100791 0.008034836
## med_int_50 0.11084350 -0.02240302 -0.09924381 -0.212842464
## high_int_100 1.00000000 -0.04673660 -0.01905821 -0.217083104
## pop_ct_1000 -0.04673660 1.00000000 -0.11493136 -0.059119958
## dist_m_cafos -0.01905821 -0.11493136 1.00000000 0.619860480
## dist_m_og_wells -0.21708310 -0.05911996 0.61986048 1.000000000
## dist_m_wwtp -0.30637142 -0.20278780 0.15352999 0.350733252
## aadt_100 0.71529092 -0.07658874 0.10155682 -0.060900696
## aadt_250 0.60158329 -0.10971668 0.23112012 0.109442838
## aadt_1000 0.40408074 0.01724558 0.12765618 0.091175296
## dist_m_wwtp aadt_100 aadt_250 aadt_1000
## tree_cover_100 0.14009235 -0.17607922 -0.007357196 0.08574999
## impervious_2500 -0.49234185 0.40911376 0.427426070 0.43050282
## low_int_100 -0.07523845 -0.25591846 -0.261080670 -0.11257695
## med_int_50 0.02725541 0.06413761 0.031971369 -0.19807436
## high_int_100 -0.30637142 0.71529092 0.601583294 0.40408074
## pop_ct_1000 -0.20278780 -0.07658874 -0.109716676 0.01724558
## dist_m_cafos 0.15352999 0.10155682 0.231120116 0.12765618
## dist_m_og_wells 0.35073325 -0.06090070 0.109442838 0.09117530
## dist_m_wwtp 1.00000000 -0.20875176 -0.024182445 -0.05838643
## aadt_100 -0.20875176 1.00000000 0.806259695 0.50962569
## aadt_250 -0.02418244 0.80625970 1.000000000 0.63521754
## aadt_1000 -0.05838643 0.50962569 0.635217540 1.00000000#' NO2, and smoke spatiotemporal predictors
#' NO2 at each site estimated using IDW
bc_st_no2 <- dplyr::select(lur_data, site_id, st_week, idw_no2) %>%
distinct() %>%
arrange(st_week) %>%
pivot_wider(id_cols = st_week,
names_from = site_id, values_from = idw_no2) %>%
as.data.frame()
rownames(bc_st_no2) <- bc_st_no2$st_week
bc_st_no2$st_week <- NULL
bc_st_no2 <- as.matrix(bc_st_no2)
no2_weeks <- rownames(bc_st_no2)
tail(no2_weeks)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"setdiff(bc_weeks, no2_weeks)## character(0)#' Smoke day indicator variable based on WFS paper method (see Brey and Fischer, 2016)
#' based on any smoke variable in the area using the 2 sd cut-off (area_smoke_2sd)
bc_st_smk <- dplyr::select(lur_data, site_id, st_week, area_smoke_2sd) %>%
distinct() %>%
arrange(st_week) %>%
pivot_wider(id_cols = st_week,
names_from = site_id, values_from = area_smoke_2sd) %>%
as.data.frame()
rownames(bc_st_smk) <- bc_st_smk$st_week
bc_st_smk$st_week <- NULL
bc_st_smk <- as.matrix(bc_st_smk)
smk_weeks <- rownames(bc_st_smk)
tail(smk_weeks)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"setdiff(bc_weeks, smk_weeks)## character(0)Create a new denver.data object.
denver.data <- createSTdata(obs = bc_obs2,
covars = bc_sp_cov,
SpatioTemporal = list(bc_st_no2 = bc_st_no2,
bc_st_smk = bc_st_smk))
D <- createDataMatrix(denver.data)
denver.data## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 231 (observed: 231)
## No. obs: 1021
##
## Only constant temporal trend,with dates:
## 2016-02-08 to 2020-11-30
##
## 18 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "high_int_100"
## [9] "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells" "dist_m_wwtp"
## [13] "aadt_100" "aadt_250" "aadt_1000" "x"
## [17] "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-062.2 Creating the NO2, BC, and PM data objects
Next we need to combine monitoring data for NO2, BC, and PM2.5 in order to calculate the time trends.
First up: NO2
Plot the time trends for NO2 at the monitors
Original units (ppb)
no2_ts <- ggplot(data = no2_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
no2_ts## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Log-transformed NO2
no2_ts2 <- ggplot(data = no2_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
no2_ts2## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data. Use the log-transformed NO2 here
#' NO2 observations
no2_obs <- no2_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(no2 = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = no2) %>%
as.data.frame() %>%
arrange(st_week)
head(no2_obs)head(no2_obs$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(no2_obs$st_week)## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"# NO2 SP object
no2_data2 <- read_csv(here::here("Data", "Monitor_NO2_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_no2")) %>%
filter(monitor_id %in% colnames(no2_obs))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.no2_coords <- do.call(rbind, st_geometry(no2_data2)) %>% as_tibble()
names(no2_coords) <- c("x", "y")
no2_sp_lonlat <- no2_data2 %>%
st_transform(crs = ll_wgs84)
no2_coords2 <- do.call(rbind, st_geometry(no2_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
no2_sp_cov <- st_set_geometry(no2_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
no2_sp_cov <- bind_cols(no2_sp_cov, no2_coords, no2_coords2) %>%
as.data.frame() %>%
distinct()Next, BC at the central site monitor
#' BC central site concentrations
bc_cent_data <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
filter(!is.na(Arithmetic_Mean)) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(County_Code = str_sub(monitor_id, start = 3, end = 5)) %>%
filter(County_Code %in% counties) %>%
arrange(Date_Local, monitor_id) %>%
mutate(Arithmetic_Mean = ifelse(Arithmetic_Mean <= 0.005, NA, Arithmetic_Mean))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Pollutant_Standard = col_logical(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_logical(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.head(bc_cent_data$Date_Local)## [1] "2016-02-08" "2016-02-09" "2016-02-10" "2016-02-11" "2016-02-12"
## [6] "2016-02-13"tail(bc_cent_data$Date_Local)## [1] "2020-11-25" "2020-11-26" "2020-11-27" "2020-11-28" "2020-11-29"
## [6] "2020-11-30"#' Add an identifier to differentiate these measurements from collocated pollutants at the same site
bc_cent_data$monitor_id <- paste0(bc_cent_data$monitor_id, "_bc")
unique(bc_cent_data$monitor_id)## [1] "080310027_bc"Plot the time trends for BC at the central site monitor
Original units (ug/m3)
bc_ts <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
#facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
bc_ts## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Log-transformed BC
bc_ts2 <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
#facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
bc_ts2## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data
Use log-transformed data here
#' BC observations
bc_cent_obs <- bc_cent_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(bc = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = bc) %>%
as.data.frame() %>%
arrange(st_week)
head(bc_cent_obs)head(bc_cent_obs$st_week)## [1] "2016-02-08" "2016-02-15" "2016-02-22" "2016-02-29" "2016-03-07"
## [6] "2016-03-14"tail(bc_cent_obs$st_week)## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"# BC Central Site SP object
bc_cent_data2 <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_bc")) %>%
filter(monitor_id %in% colnames(bc_cent_obs))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Pollutant_Standard = col_logical(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_logical(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.bc_cent_coords <- do.call(rbind, st_geometry(bc_cent_data2)) %>% as_tibble()
names(bc_cent_coords) <- c("x", "y")
bc_cent_sp_lonlat <- bc_cent_data2 %>%
st_transform(crs = ll_wgs84)
bc_cent_coords2 <- do.call(rbind, st_geometry(bc_cent_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
bc_cent_sp_cov <- st_set_geometry(bc_cent_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
bc_cent_sp_cov <- bind_cols(bc_cent_sp_cov, bc_cent_coords, bc_cent_coords2) %>%
as.data.frame() %>%
distinct()Lastly, PM2.5 from select monitors with long-term records
Plot the time trends for PM2.5 at the monitors
Original units (ug/m3)
pm_ts <- ggplot(data = pm_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
pm_ts## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Log-transformed PM2.5
pm_ts2 <- ggplot(data = pm_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
pm_ts2## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data
Use the log-transformed data here
#' PM2.5 observations
pm_obs <- pm_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(pm = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = pm) %>%
as.data.frame() %>%
arrange(st_week)
head(pm_obs)head(pm_obs$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(pm_obs$st_week)## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"# PM SP object
pm_data2 <- read_csv(here::here("Data", "Monitor_PM_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_pm")) %>%
filter(monitor_id %in% colnames(pm_obs))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_double(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.pm_coords <- do.call(rbind, st_geometry(pm_data2)) %>% as_tibble()
names(pm_coords) <- c("x", "y")
pm_sp_lonlat <- pm_data2 %>%
st_transform(crs = ll_wgs84)
pm_coords2 <- do.call(rbind, st_geometry(pm_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
pm_sp_cov <- st_set_geometry(pm_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
pm_sp_cov <- bind_cols(pm_sp_cov, pm_coords, pm_coords2) %>%
as.data.frame() %>%
distinct()Combine the observations and the spatial covariates Scale the pollutant measurements
#' Join log-transformed observations
pol_obs <- left_join(no2_obs, bc_cent_obs, by = "st_week") %>%
left_join(pm_obs, by = "st_week")
glimpse(pol_obs)## Rows: 622
## Columns: 18
## $ st_week <chr> "2008-12-29", "2009-01-05", "2009-01-12", "2009-01-19…
## $ `080013001_no2` <dbl> 2.4313422, 2.3337153, 3.0953820, 2.4960668, 2.6445846…
## $ `080310002_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310026_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_bc` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080010006_pm` <dbl> 1.394118, 1.247789, 1.934844, 2.298456, 1.598321, 1.7…
## $ `080010008_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080050005_pm` <dbl> 1.2490759, 0.6817687, 1.3697744, 1.7730816, 1.1939225…
## $ `080130003_pm` <dbl> 0.9679299, 0.8243293, 2.1047287, 1.7212142, 1.9059908…
## $ `080130012_pm` <dbl> 1.0919008, 0.8908546, 1.2511276, 1.2703657, 0.8047190…
## $ `080310002_pm` <dbl> 1.335209, 1.188381, 1.767997, 2.027844, 1.486213, 1.7…
## $ `080310023_pm` <dbl> 1.661197, 1.513911, 1.919262, 1.996528, 1.909773, 1.7…
## $ `080310025_pm` <dbl> 1.1537863, 0.9359011, 1.3262687, 1.9969858, 1.6521301…
## $ `080310026_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…rownames(pol_obs) <- pol_obs$st_week
pol_obs$st_week <- NULL
pol_obs <- as.matrix(pol_obs)
#' Scale the measurements to avoid issues where units are different
#' Remember! These have been log-transformed before scaling
#' Scale the data by columns (monitors)
rnames <- rownames(pol_obs)
cnames <- colnames(pol_obs)
pol_obs <- apply(pol_obs, 2, scale)
rownames(pol_obs) <- rnames
colnames(pol_obs) <- cnames
head(rownames(pol_obs))## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(rownames(pol_obs))## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"head(colnames(pol_obs))## [1] "080013001_no2" "080310002_no2" "080310026_no2" "080310027_no2"
## [5] "080310028_no2" "080310027_bc"class(pol_obs)## [1] "matrix" "array"dim(pol_obs)## [1] 622 17#' SP variables
pol_sp_cov <- bind_rows(no2_sp_cov, bc_cent_sp_cov, pm_sp_cov)
pol_sp_covCreate the data object
pol_STdata <- createSTdata(obs = pol_obs,
covars = pol_sp_cov)
print(pol_STdata)## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Only constant temporal trend,with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.2.3 Generate the temporal trend
The cross validation results for generating the time trends suggest 2 basis functions is the way to go, but let’s see if that holds up with our data
D_pol <- createDataMatrix(pol_STdata)
#D_pol
n_years <- length(2009:2020)
n_years*4## [1] 48n_years*8## [1] 96#' Trying 4 df per year
pol.SVD.cv.4py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 4)
print(pol.SVD.cv.4py)## Result of SVDsmoothCV, average of CV-statistics:
## MSE R2 AIC BIC
## n.basis.0 0.9975610 0.0000000 0.9987781 5.095695
## n.basis.1 0.8013682 0.1967023 -96.6124965 -88.418663
## n.basis.2 0.6358644 0.3625750 -225.3897791 -213.099029
## n.basis.3 0.6297167 0.3687387 -229.4297802 -213.042114
## n.basis.4 0.6084604 0.3900698 -240.5027300 -220.018147plot(pol.SVD.cv.4py)
#' Trying 8 df per year
pol.SVD.cv.8py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 8)
print(pol.SVD.cv.8py)## Result of SVDsmoothCV, average of CV-statistics:
## MSE R2 AIC BIC
## n.basis.0 0.9975610 0.0000000 0.9987781 5.095695
## n.basis.1 0.7630671 0.2350911 -120.2211615 -112.027328
## n.basis.2 0.5825979 0.4159654 -265.3105795 -253.019830
## n.basis.3 0.5739736 0.4246129 -271.3027690 -254.915103
## n.basis.4 0.5615078 0.4371173 -278.8000243 -258.315442plot(pol.SVD.cv.8py)
2.4 Plot the time trends
I tried using 1 and 2 basis functions for the time trends. I also checked out 4 and 8 df per year
2.4.1 1 basis function
First, see what the time trends for the monitoring data look like with 1 basis function.
4 df per year
pol_STdata1 <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 4)## Replacing existing trend.pol_STdata1## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata1$trend)tail(pol_STdata1$trend)plot(pol_STdata1$trend$date, pol_STdata1$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1$trend$date, pol_STdata1$atrend$V1, col = 1)
8 df per year
pol_STdata1a <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 8)## Replacing existing trend.pol_STdata1a## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata1a$trend)tail(pol_STdata1a$trend)plot(pol_STdata1a$trend$date, pol_STdata1a$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1a$trend$date, pol_STdata1a$atrend$V1, col = 1)
2.4.2 2 basis functions
Next, let’s see what 2 basis functions look like:
4 df per year
pol_STdata2 <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 4)## Replacing existing trend.pol_STdata2## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata2$trend)tail(pol_STdata2$trend)plot(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2)
8 df per year
pol_STdata2a <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 8)## Replacing existing trend.pol_STdata2a## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata2a$trend)tail(pol_STdata2a$trend)plot(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2)
2.5 Update the Denver data set
The next step is to update the trend data for the denver.data object and see if these trend data fit our observations.
Comparing the plots using one basis function vs. two basis functions, it looks like using two basis functions might be better than 1. There was some residual autocorrelation in the central site monitor data when using just one basis function (see the fourth plot in the first series of plots below). Using 8 degrees of freedom did not help.
Using two basis functions helped somewhat to reduce this autocorrelation, but it still persists. There are still some noticable trends in the residuals for the central monitoring site
2.5.1 1 basis function, 4 df per year
First update the trend data, then plot the trends for three distributed sites and the central BC monitoring site.
denver.data2.1 <- denver.data
denver.data2.1$trend <- pol_STdata1$trend
print(denver.data2.1)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 18 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "high_int_100"
## [9] "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells" "dist_m_wwtp"
## [13] "aadt_100" "aadt_250" "aadt_1000" "x"
## [17] "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.1, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_1")
plot(denver.data2.1, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_10")
plot(denver.data2.1, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_43")
plot(denver.data2.1, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.1, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_53")
plot(denver.data2.1, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.1, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="central")
plot(denver.data2.1, "pacf", ID="central")
2.5.2 1 basis function, 8 df per year
denver.data2.1a <- denver.data
denver.data2.1a$trend <- pol_STdata1a$trend
print(denver.data2.1a)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 18 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "high_int_100"
## [9] "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells" "dist_m_wwtp"
## [13] "aadt_100" "aadt_250" "aadt_1000" "x"
## [17] "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.1a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_1")
plot(denver.data2.1a, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_10")
plot(denver.data2.1a, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_43")
plot(denver.data2.1a, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.1a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_53")
plot(denver.data2.1a, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.1a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="central")
plot(denver.data2.1a, "pacf", ID="central")
2.5.3 2 basis functions, 4 df per year
denver.data2.2 <- denver.data
denver.data2.2$trend <- pol_STdata2$trend
print(denver.data2.2)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 18 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "high_int_100"
## [9] "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells" "dist_m_wwtp"
## [13] "aadt_100" "aadt_250" "aadt_1000" "x"
## [17] "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.2, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_1")
plot(denver.data2.2, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_10")
plot(denver.data2.2, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_10")
plot(denver.data2.2, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.2, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_43")
plot(denver.data2.2, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.2, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_53")
plot(denver.data2.2, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.2, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="central")
plot(denver.data2.2, "pacf", ID="central")
2.5.4 2 basis functions, 8 df per year
denver.data2.2a <- denver.data
denver.data2.2a$trend <- pol_STdata2a$trend
print(denver.data2.2a)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 18 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "high_int_100"
## [9] "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells" "dist_m_wwtp"
## [13] "aadt_100" "aadt_250" "aadt_1000" "x"
## [17] "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.2a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_1")
plot(denver.data2.2a, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_10")
plot(denver.data2.2a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_10")
plot(denver.data2.2a, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.2a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_43")
plot(denver.data2.2a, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.2a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_53")
plot(denver.data2.2a, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.2a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="central")
plot(denver.data2.2a, "pacf", ID="central")
3 Testing possible models
Building on the results of the previous modeling efforts, I’m going to update “Model 4.2”, which used an ’exp covariance structure for the error term and iid for the beta fields.
3.1 Model A: 1 temporal trend (basis) function (4 df per year)
3.1.1 Create the model object
For this version of the model, use iid for cov.beta (beta, beta1) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.
denver.data.A <- denver.data2.1
names(denver.data.A$covars)## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "high_int_100"
## [9] "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells" "dist_m_wwtp"
## [13] "aadt_100" "aadt_250" "aadt_1000" "x"
## [17] "y" "type"LUR.A <- list(covar_fun, covar_fun)
cov.beta.A <- list(covf = c("iid", "iid"), nugget = T)
cov.nu.A <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.A <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")
denver.model.A <- createSTmodel(denver.data.A, LUR = LUR.A,
ST = c("bc_st_no2", "bc_st_smk"),
cov.beta = cov.beta.A, cov.nu = cov.nu.A,
locations = locations.A)## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.denver.model.A## STmodel-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## Models for the beta-fields are:
## $const
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## high_int_100 + pop_ct_1000 + dist_m_cafos + dist_m_og_wells +
## dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## $V1
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## high_int_100 + pop_ct_1000 + dist_m_cafos + dist_m_og_wells +
## dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## Covariance model for the beta-field(s):
## Covariance type(s): iid, iid
## Nugget: Yes, Yes
## Covariance model for the nu-field(s):
## Covariance type: exp
## Nugget: ~1
## Random effect: No
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-063.1.2 Estimate model parameters
- nugget values can range from -1 to -6
- range values can range from 0 to 4
- try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
- make sure the starting values are pretty different
set.seed(1000)
names <- loglikeSTnames(denver.model.A, all=FALSE)
names## [1] "log.nugget.const.iid" "log.nugget.V1.iid"
## [3] "nu.log.range.exp" "nu.log.sill.exp"
## [5] "nu.log.nugget.(Intercept).exp"# x.init.A <- cbind(c(0, 0, 0, 0, 0),
# c(-1, -1, 0, -1, -1),
# c(-1, -1, 0, -5, -1),
# c(-5, -5, 0, -1, -5),
# c(-5, -5, 0, -5, -5),
# c(-1, -1, 2, -1, -1),
# c(-1, -1, 2, -5, -1),
# c(-5, -5, 2, -1, -5),
# c(-5, -5, 2, -5, -5),
# c(-1, -1, 4, -1, -1),
# c(-1, -1, 4, -5, -1),
# c(-5, -5, 4, -1, -5),
# c(-5, -5, 4, -5, -5))
x.init.A <- cbind(c(0, 0, 0, 0, 0),
c(-1, -1, 4, -5, -1),
c(-5, -5, 4, -1, -5),
c(-1, -1, 8, -5, -1),
c(-5, -5, 8, -1, -5),
c(-1, -1, 12, -5, -1),
c(-5, -5, 12, -1, -5),
c(-7, -7, 12, -3, -3))
rownames(x.init.A) <- loglikeSTnames(denver.model.A, all=FALSE)
x.init.A## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## log.nugget.const.iid 0 -1 -5 -1 -5 -1 -5 -7
## log.nugget.V1.iid 0 -1 -5 -1 -5 -1 -5 -7
## nu.log.range.exp 0 4 4 8 8 12 12 12
## nu.log.sill.exp 0 -5 -1 -5 -1 -5 -1 -3
## nu.log.nugget.(Intercept).exp 0 -1 -5 -1 -5 -1 -5 -3#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.A <- estimate.STmodel(denver.model.A, x.init.A)## Optimisation using starting value 1/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= 459.08 |proj g|= 15
## At iterate 10 f = -1308.6 |proj g|= 0.87075
##
## iterations 19
## function evaluations 32
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00264401
## final function value -1308.8
##
## F = -1308.8
## final value -1308.797217
## converged## Optimisation using starting value 2/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -331.86 |proj g|= 14
## At iterate 10 f = -1305 |proj g|= 1.8915
## ys=-2.226e-01 -gs= 3.952e-02, BFGS update SKIPPED
## At iterate 20 f = -1306.1 |proj g|= 4.6671
## At iterate 30 f = -1308.8 |proj g|= 0.22097
## At iterate 40 f = -1433.2 |proj g|= 20.29
## At iterate 50 f = -1570.9 |proj g|= 4.1779
##
## iterations 54
## function evaluations 85
## segments explored during Cauchy searches 57
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0297467
## final function value -1570.92
##
## F = -1570.92
## final value -1570.918992
## converged## Optimisation using starting value 3/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -458.59 |proj g|= 14
## At iterate 10 f = -1303.6 |proj g|= 1.8774
## At iterate 20 f = -1308.5 |proj g|= 5.5309
## At iterate 30 f = -1570.9 |proj g|= 2.5698
##
## iterations 34
## function evaluations 54
## segments explored during Cauchy searches 35
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.112304
## final function value -1570.92
##
## F = -1570.92
## final value -1570.920464
## converged## Optimisation using starting value 4/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -332.17 |proj g|= 14
## At iterate 10 f = -1568.5 |proj g|= 7.865
## ys=-2.281e-01 -gs= 6.765e-01, BFGS update SKIPPED
## At iterate 20 f = -1575.1 |proj g|= 0.49892
## At iterate 30 f = -1575.2 |proj g|= 9.8017e-05
##
## iterations 30
## function evaluations 51
## segments explored during Cauchy searches 33
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 9.80167e-05
## final function value -1575.2
##
## F = -1575.2
## final value -1575.197734
## converged## Optimisation using starting value 5/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -546.64 |proj g|= 14
## At iterate 10 f = -1575 |proj g|= 0.99606
## At iterate 20 f = -1575.2 |proj g|= 0.022396
##
## iterations 24
## function evaluations 44
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0100915
## final function value -1575.2
##
## F = -1575.2
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1575.197751
## converged## Optimisation using starting value 6/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -335.65 |proj g|= 14
## At iterate 10 f = -1575 |proj g|= 6.9895
##
## iterations 18
## function evaluations 42
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0331148
## final function value -1575.2
##
## F = -1575.2
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1575.197763
## converged## Optimisation using starting value 7/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1296.6 |proj g|= 14
## At iterate 10 f = -1575.2 |proj g|= 0.30063
## At iterate 20 f = -1575.2 |proj g|= 0.044572
##
## iterations 25
## function evaluations 43
## segments explored during Cauchy searches 29
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00497526
## final function value -1575.2
##
## F = -1575.2
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1575.197749
## converged## Optimisation using starting value 8/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1228.6 |proj g|= 12
## At iterate 10 f = -1575.2 |proj g|= 1.5217
## At iterate 20 f = -1575.2 |proj g|= 0.11938
##
## iterations 21
## function evaluations 27
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.100619
## final function value -1575.2
##
## F = -1575.2
## final value -1575.197742
## convergedprint(est.denver.model.A)## Optimisation for STmodel with 8 starting points.
## Results: 5 converged, 3 not converged, 0 failed.
## Best result for starting point 6, optimisation has converged
##
## No fixed parameters.
##
## Estimated parameters for all starting point(s):
## [,1] [,2] [,3]
## gamma.bc_st_no2 0.0188524702 0.0122849478 0.0122935795
## gamma.bc_st_smk -0.0463102718 -0.0241238841 -0.0241068538
## alpha.const.(Intercept) -0.1467826748 0.0397887383 0.0396047027
## alpha.const.tree_cover_100 -0.0005338163 0.0045194887 0.0045150989
## alpha.const.impervious_2500 0.0402633123 0.0269217938 0.0269125398
## alpha.const.low_int_100 -0.0041473553 0.0078694534 0.0078669871
## alpha.const.med_int_50 -0.0030576675 0.0034919877 0.0034934544
## alpha.const.high_int_100 0.0104424347 0.0167840080 0.0167813996
## alpha.const.pop_ct_1000 -0.0133634318 -0.0154765220 -0.0154708819
## alpha.const.dist_m_cafos -0.0426982983 -0.0345475238 -0.0345532847
## alpha.const.dist_m_og_wells 0.0105183299 0.0044043124 0.0044106846
## alpha.const.dist_m_wwtp 0.0027145445 0.0046942753 0.0046913105
## alpha.const.aadt_100 0.0358726241 0.0345600491 0.0345600984
## alpha.const.aadt_250 0.0081137884 0.0053358686 0.0053291095
## alpha.const.aadt_1000 -0.0013442311 0.0010965896 0.0010994992
## alpha.V1.(Intercept) 0.1859389126 0.2191504409 0.2191121360
## alpha.V1.tree_cover_100 -0.0192119023 -0.0149304598 -0.0149367227
## alpha.V1.impervious_2500 0.0054979907 -0.0122830864 -0.0122900148
## alpha.V1.low_int_100 -0.0094298447 -0.0017948749 -0.0017992738
## alpha.V1.med_int_50 -0.0238386525 -0.0185573635 -0.0185599707
## alpha.V1.high_int_100 0.0084869198 0.0131483199 0.0131370621
## alpha.V1.pop_ct_1000 -0.0030220235 -0.0011844435 -0.0011858298
## alpha.V1.dist_m_cafos 0.0046301267 0.0008554407 0.0008447220
## alpha.V1.dist_m_og_wells -0.0052729031 0.0003980968 0.0003945947
## alpha.V1.dist_m_wwtp -0.0070082982 -0.0090806239 -0.0090877200
## alpha.V1.aadt_100 -0.0470195439 -0.0429313086 -0.0429215815
## alpha.V1.aadt_250 0.0094461581 0.0072705831 0.0072748879
## alpha.V1.aadt_1000 -0.0075137095 -0.0052008772 -0.0051994064
## log.nugget.const.iid -14.9390328408 -14.9992686779 -14.9996048867
## log.nugget.V1.iid -15.0000000000 -14.9998091565 -14.4615169601
## nu.log.range.exp 0.0000000000 12.4788563301 12.4746448599
## nu.log.sill.exp -4.3851285631 -3.0229529932 -3.0235904960
## nu.log.nugget.(Intercept).exp -4.1355654291 -4.7649855334 -4.7650905133
## [,4] [,5] [,6]
## gamma.bc_st_no2 0.0123558241 0.0123543528 0.0123549645
## gamma.bc_st_smk -0.0232712280 -0.0232745349 -0.0232715834
## alpha.const.(Intercept) 0.0336966099 0.0337313445 0.0337142820
## alpha.const.tree_cover_100 0.0036297032 0.0036316171 0.0036310080
## alpha.const.impervious_2500 0.0264100677 0.0264128814 0.0264127240
## alpha.const.low_int_100 0.0068722755 0.0068739596 0.0068731018
## alpha.const.med_int_50 0.0023638264 0.0023642410 0.0023634899
## alpha.const.high_int_100 0.0148465408 0.0148493142 0.0148476278
## alpha.const.pop_ct_1000 -0.0141987916 -0.0142012896 -0.0142000047
## alpha.const.dist_m_cafos -0.0345176581 -0.0345169534 -0.0345169908
## alpha.const.dist_m_og_wells 0.0042218048 0.0042209047 0.0042207274
## alpha.const.dist_m_wwtp 0.0038971481 0.0038985755 0.0038979913
## alpha.const.aadt_100 0.0342280959 0.0342283265 0.0342283152
## alpha.const.aadt_250 0.0056562176 0.0056572994 0.0056576185
## alpha.const.aadt_1000 0.0016665371 0.0016648768 0.0016654639
## alpha.V1.(Intercept) 0.2230533103 0.2230646044 0.2230648289
## alpha.V1.tree_cover_100 -0.0152000058 -0.0151981517 -0.0151981960
## alpha.V1.impervious_2500 -0.0122587079 -0.0122555773 -0.0122553716
## alpha.V1.low_int_100 -0.0030485913 -0.0030462096 -0.0030468825
## alpha.V1.med_int_50 -0.0171043970 -0.0171052150 -0.0171043231
## alpha.V1.high_int_100 0.0125897783 0.0125931796 0.0125933398
## alpha.V1.pop_ct_1000 -0.0011045139 -0.0011050087 -0.0011048460
## alpha.V1.dist_m_cafos 0.0005137008 0.0005164112 0.0005167916
## alpha.V1.dist_m_og_wells -0.0001074939 -0.0001061171 -0.0001065371
## alpha.V1.dist_m_wwtp -0.0074615204 -0.0074618741 -0.0074602339
## alpha.V1.aadt_100 -0.0404261690 -0.0404304843 -0.0404299202
## alpha.V1.aadt_250 0.0073599070 0.0073580353 0.0073585462
## alpha.V1.aadt_1000 -0.0049644872 -0.0049652701 -0.0049656871
## log.nugget.const.iid -7.7887866529 -7.7903228850 -7.7891553958
## log.nugget.V1.iid -7.3764905741 -7.3781315180 -7.3770135459
## nu.log.range.exp 12.4461656570 12.4469412808 12.4468659837
## nu.log.sill.exp -3.0077198203 -3.0077604407 -3.0078350936
## nu.log.nugget.(Intercept).exp -4.8617409014 -4.8615308790 -4.8614493842
## [,7] [,8]
## gamma.bc_st_no2 0.0123550170 0.0123524584
## gamma.bc_st_smk -0.0232728166 -0.0232814448
## alpha.const.(Intercept) 0.0337147425 0.0337747113
## alpha.const.tree_cover_100 0.0036305136 0.0036279582
## alpha.const.impervious_2500 0.0264114906 0.0264106474
## alpha.const.low_int_100 0.0068728836 0.0068710833
## alpha.const.med_int_50 0.0023637901 0.0023631000
## alpha.const.high_int_100 0.0148474420 0.0148445657
## alpha.const.pop_ct_1000 -0.0141997619 -0.0141983798
## alpha.const.dist_m_cafos -0.0345171654 -0.0345161465
## alpha.const.dist_m_og_wells 0.0042212001 0.0042220904
## alpha.const.dist_m_wwtp 0.0038977249 0.0038964101
## alpha.const.aadt_100 0.0342282801 0.0342296478
## alpha.const.aadt_250 0.0056569128 0.0056563599
## alpha.const.aadt_1000 0.0016658626 0.0016656200
## alpha.V1.(Intercept) 0.2230598904 0.2230632243
## alpha.V1.tree_cover_100 -0.0151990340 -0.0152006575
## alpha.V1.impervious_2500 -0.0122570943 -0.0122628563
## alpha.V1.low_int_100 -0.0030475585 -0.0030509040
## alpha.V1.med_int_50 -0.0171044875 -0.0171038844
## alpha.V1.high_int_100 0.0125916527 0.0125904615
## alpha.V1.pop_ct_1000 -0.0011046799 -0.0011019352
## alpha.V1.dist_m_cafos 0.0005151901 0.0005172158
## alpha.V1.dist_m_og_wells -0.0001069144 -0.0001118214
## alpha.V1.dist_m_wwtp -0.0074611425 -0.0074582043
## alpha.V1.aadt_100 -0.0404281327 -0.0404197671
## alpha.V1.aadt_250 0.0073590565 0.0073591319
## alpha.V1.aadt_1000 -0.0049649754 -0.0049643332
## log.nugget.const.iid -7.7891990039 -7.7874513701
## log.nugget.V1.iid -7.3769712527 -7.3729570707
## nu.log.range.exp 12.4466055842 12.4471890786
## nu.log.sill.exp -3.0077180036 -3.0066509945
## nu.log.nugget.(Intercept).exp -4.8616295839 -4.8617088663
##
## Function value(s):
## [1] 1308.797 1570.919 1570.920 1575.198 1575.198 1575.198 1575.198 1575.1983.1.3 Cross-validation
Define the CV groups
set.seed(1000)
site_idsA <- colnames(denver.model.A$D.beta)[which(str_detect(colnames(denver.model.A$D.beta), "d_") == T)]
unique(colnames(bc_obs2))## [1] "central" "d_1" "d_10" "d_11" "d_12" "d_13" "d_14"
## [8] "d_15" "d_16" "d_17" "d_18" "d_19" "d_2" "d_20"
## [15] "d_21" "d_22" "d_23" "d_24" "d_25" "d_26" "d_27"
## [22] "d_28" "d_29" "d_3" "d_30" "d_31" "d_32" "d_33"
## [29] "d_34" "d_35" "d_36" "d_37" "d_38" "d_39" "d_4"
## [36] "d_40" "d_41" "d_42" "d_43" "d_44" "d_45" "d_46"
## [43] "d_47" "d_48" "d_49" "d_5" "d_50" "d_51" "d_52"
## [50] "d_53" "d_54" "d_55" "d_56" "d_57" "d_58" "d_59"
## [57] "d_6" "d_60" "d_61" "d_62" "d_63" "d_64" "d_65"
## [64] "d_66" "d_67" "d_68" "d_7" "d_8" "d_9"Ind.cv.A <- createCV(denver.model.A, groups = 10, min.dist = .1,
subset = site_idsA)
ID.cv.A <- sapply(split(denver.model.A$obs$ID, Ind.cv.A), unique)
print(sapply(ID.cv.A, length))## 0 1 2 3 4 5 6 7 8 9 10
## 1 7 7 7 6 6 6 6 6 6 6table(Ind.cv.A)## Ind.cv.A
## 0 1 2 3 4 5 6 7 8 9 10
## 231 75 80 97 73 64 93 80 81 66 81I.col.A <- apply(sapply(ID.cv.A, function(x) denver.model.A$locations$ID%in% x), 1,
function(x) if(sum(x)==1) which(x) else 0)
names(I.col.A) <- denver.model.A$locations$ID
print(I.col.A)## central d_1 d_10 d_11 d_12 d_13 d_14 d_15 d_16 d_17
## 1 6 2 3 4 4 4 9 4 10
## d_18 d_19 d_2 d_20 d_21 d_22 d_23 d_25 d_26 d_28
## 8 5 11 3 2 11 4 9 11 9
## d_29 d_3 d_31 d_32 d_33 d_34 d_35 d_36 d_37 d_38
## 2 4 6 8 5 5 10 8 3 8
## d_39 d_4 d_40 d_41 d_42 d_43 d_44 d_45 d_46 d_47
## 6 7 7 3 7 8 11 3 9 5
## d_48 d_49 d_5 d_50 d_51 d_52 d_53 d_54 d_55 d_56
## 9 4 6 9 10 3 7 5 7 6
## d_57 d_58 d_6 d_60 d_61 d_62 d_63 d_64 d_65 d_66
## 2 2 11 5 2 11 10 3 8 6
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 7 10 2 10 0 0 0 0 0par(mfrow=c(1,1))
plot(denver.model.A$locations$long,
denver.model.A$locations$lat,
pch=23+floor(I.col.A/max(I.col.A)+.5), bg=I.col.A,
xlab="Longitude", ylab="Latitude")
map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)## [[1]]
## NULLID starting conditions using previous model without CV:
x.init.A.cv <- coef(est.denver.model.A, pars="cov")[,c("par","init")]
x.init.A.cvRun the model with cross validation
est.denver.A.cv <- estimateCV(denver.model.A, x.init.A.cv, Ind.cv.A)##
## ***************************## Estimation of CV-set 1/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1431.1 |proj g|= 10.81
## At iterate 10 f = -1432.1 |proj g|= 0.21349
## At iterate 20 f = -1432.1 |proj g|= 0.26792
##
## iterations 24
## function evaluations 33
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0427763
## final function value -1432.08
##
## F = -1432.08
## final value -1432.078392
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -312.21 |proj g|= 14
## At iterate 10 f = -1431.9 |proj g|= 5.3771
## At iterate 20 f = -1432.1 |proj g|= 0.4494
##
## iterations 25
## function evaluations 35
## segments explored during Cauchy searches 28
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0398318
## final function value -1432.08
##
## F = -1432.08
## final value -1432.078408
## converged## ##
## ***************************## Estimation of CV-set 2/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1434 |proj g|= 1.401
## At iterate 10 f = -1434.1 |proj g|= 0.079292
## At iterate 20 f = -1434.1 |proj g|= 0.0092636
##
## iterations 20
## function evaluations 36
## segments explored during Cauchy searches 20
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00926357
## final function value -1434.14
##
## F = -1434.14
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1434.139854
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -312.02 |proj g|= 14
## At iterate 10 f = -1434.1 |proj g|= 0.039452
##
## iterations 15
## function evaluations 30
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0645088
## final function value -1434.14
##
## F = -1434.14
## final value -1434.139857
## converged## ##
## ***************************## Estimation of CV-set 3/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1420.1 |proj g|= 10.139
## At iterate 10 f = -1420.9 |proj g|= 0.39574
##
## iterations 12
## function evaluations 21
## segments explored during Cauchy searches 14
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0802538
## final function value -1420.91
##
## F = -1420.91
## final value -1420.907659
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -304.38 |proj g|= 14
## At iterate 10 f = -1420.1 |proj g|= 10.114
##
## iterations 19
## function evaluations 37
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0381455
## final function value -1420.91
##
## F = -1420.91
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1420.907641
## converged## ##
## ***************************## Estimation of CV-set 4/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1424.4 |proj g|= 17.069
## At iterate 10 f = -1425.3 |proj g|= 0.15949
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 17
## function evaluations 42
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00646122
## final function value -1425.26
##
## F = -1425.26
## final value -1425.262876
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -311.89 |proj g|= 14
## At iterate 10 f = -1425 |proj g|= 3.896
## At iterate 20 f = -1425.3 |proj g|= 0.13831
##
## iterations 24
## function evaluations 43
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00885772
## final function value -1425.26
##
## F = -1425.26
## final value -1425.262862
## converged## ##
## ***************************## Estimation of CV-set 5/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1451.2 |proj g|= 9.2005
## At iterate 10 f = -1451.6 |proj g|= 0.79403
##
## iterations 19
## function evaluations 27
## segments explored during Cauchy searches 19
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0186199
## final function value -1451.65
##
## F = -1451.65
## final value -1451.648175
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -316.56 |proj g|= 14
## At iterate 10 f = -1451.6 |proj g|= 0.12505
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1451.648203
## stopped after 19 iterations## ##
## ***************************## Estimation of CV-set 6/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1421.6 |proj g|= 10.093
## At iterate 10 f = -1422 |proj g|= 0.24382
##
## iterations 18
## function evaluations 33
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0392072
## final function value -1421.97
##
## F = -1421.97
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1421.969444
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -305.34 |proj g|= 14
## At iterate 10 f = -1421.8 |proj g|= 1.2766
## At iterate 20 f = -1422 |proj g|= 0.039819
##
## iterations 22
## function evaluations 45
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0482915
## final function value -1421.97
##
## F = -1421.97
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1421.969411
## converged## ##
## ***************************## Estimation of CV-set 7/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1431.4 |proj g|= 7.6155
## At iterate 10 f = -1433.5 |proj g|= 3.6513
## At iterate 20 f = -1433.7 |proj g|= 0.81046
##
## iterations 28
## function evaluations 36
## segments explored during Cauchy searches 28
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0572387
## final function value -1433.73
##
## F = -1433.73
## final value -1433.728979
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -310.57 |proj g|= 14
## At iterate 10 f = -1433.3 |proj g|= 5.7091
## At iterate 20 f = -1433.7 |proj g|= 0.13641
## At iterate 30 f = -1433.7 |proj g|= 0.009134
##
## iterations 32
## function evaluations 48
## segments explored during Cauchy searches 35
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00857919
## final function value -1433.73
##
## F = -1433.73
## final value -1433.728950
## converged## ##
## ***************************## Estimation of CV-set 8/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1421.7 |proj g|= 8.2282
## At iterate 10 f = -1422.2 |proj g|= 0.97761
##
## iterations 19
## function evaluations 23
## segments explored during Cauchy searches 19
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00683795
## final function value -1422.22
##
## F = -1422.22
## final value -1422.223586
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -310.64 |proj g|= 14
## At iterate 10 f = -1421.9 |proj g|= 9.2206
## At iterate 20 f = -1422.2 |proj g|= 0.032827
##
## iterations 20
## function evaluations 40
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0328272
## final function value -1422.22
##
## F = -1422.22
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1422.223598
## converged## ##
## ***************************## Estimation of CV-set 9/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1443.4 |proj g|= 12.678
## At iterate 10 f = -1444.2 |proj g|= 0.20454
## At iterate 20 f = -1444.2 |proj g|= 0.010781
##
## iterations 22
## function evaluations 32
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.046591
## final function value -1444.18
##
## F = -1444.18
## final value -1444.177944
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -315.18 |proj g|= 14
## At iterate 10 f = -1444.2 |proj g|= 0.42178
## At iterate 20 f = -1444.2 |proj g|= 0.013834
##
## iterations 23
## function evaluations 40
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0063626
## final function value -1444.18
##
## F = -1444.18
## final value -1444.177925
## converged## ##
## ***************************## Estimation of CV-set 10/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1447.9 |proj g|= 10.139
## At iterate 10 f = -1450.2 |proj g|= 0.071665
##
## iterations 16
## function evaluations 20
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00478007
## final function value -1450.16
##
## F = -1450.16
## final value -1450.162527
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -309.55 |proj g|= 14
## At iterate 10 f = -1449.8 |proj g|= 2.9604
## At iterate 20 f = -1450.2 |proj g|= 0.41207
##
## iterations 23
## function evaluations 43
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0284068
## final function value -1450.16
##
## F = -1450.16
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1450.162578
## converged## print(est.denver.A.cv)## Cross-validation parameter estimation for STmodel
## with 10 CV-groups and 2 starting points.
## Results: 10 converged, 0 not converged.
##
## No fixed parameters.
##
## Estimated function values and convergence info:
## value convergence conv eigen.min eigen.all.min
## 1 1432.078 TRUE TRUE 0.32931870 NA
## 2 1434.140 TRUE TRUE 2.22583856 NA
## 3 1420.908 TRUE TRUE 2.59289894 NA
## 4 1425.263 TRUE TRUE 0.84473934 NA
## 5 1451.648 TRUE TRUE 0.98039663 NA
## 6 1421.969 TRUE TRUE 2.04327199 NA
## 7 1433.729 TRUE TRUE 0.05845734 NA
## 8 1422.224 TRUE TRUE 0.62572871 NA
## 9 1444.178 TRUE TRUE 1.15115963 NA
## 10 1450.163 TRUE TRUE 3.06878920 NApar(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.A, pars="all"),
plotCI((1:length(par))+.3, par, uiw=1.96*sd,
col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.A.cv, "all", boxwex=.4, col="grey", add=TRUE)
#' Save the results as an .rdata object
save(denver.data.A, denver.model.A, est.denver.model.A, est.denver.A.cv,
file = here::here("Results", "Denver_ST_Model_A.rdata"))3.1.4 Prediction using the CV model
Making predictions using the CV model. Printing out the CV summary statistics as well
pred.A.cv <- predictCV(denver.model.A, est.denver.A.cv, LTA = T)
pred.A.cv.log <- predictCV(denver.model.A, est.denver.A.cv,
LTA = T, transform="unbiased")
head(pred.A.cv$pred.all$EX)## central d_1 d_10 d_11 d_12
## 2008-12-29 NA -0.37207606 -0.234784391 -0.1795091567 -0.233196632
## 2009-01-05 NA -0.29145083 -0.163169432 -0.1102228971 -0.162492338
## 2009-01-12 NA -0.13511050 -0.017201930 0.0342488286 -0.029315156
## 2009-01-19 NA -0.15907650 -0.047677254 0.0005574218 -0.045251704
## 2009-01-26 NA -0.09824328 0.005436241 0.0515219618 0.008605429
## 2009-02-02 NA -0.01561741 0.080516360 0.1249336890 0.080130733
## d_13 d_14 d_15 d_16 d_17
## 2008-12-29 -0.26576323 -0.44241186 -0.41770498 0.2197111 -0.41489784
## 2009-01-05 -0.19860989 -0.36322039 -0.33145257 0.2645376 -0.32459735
## 2009-01-12 -0.06891953 -0.22170929 -0.17082573 0.3723043 -0.16355599
## 2009-01-19 -0.08819938 -0.22965497 -0.18788999 0.3320032 -0.17220088
## 2009-01-26 -0.03746260 -0.16833983 -0.12195115 0.3631205 -0.10224277
## 2009-02-02 0.03123185 -0.09004848 -0.03513389 0.4140159 -0.01304501
## d_18 d_19 d_2 d_20 d_21
## 2008-12-29 -0.30621502 -0.277865155 0.2949866 -0.27593821 -0.31282340
## 2009-01-05 -0.23035094 -0.198319634 0.3418046 -0.20014273 -0.23644923
## 2009-01-12 -0.07896986 -0.040553719 0.4646358 -0.04927932 -0.08580847
## 2009-01-19 -0.10705093 -0.068753445 0.4092439 -0.07684215 -0.11180290
## 2009-01-26 -0.05036466 -0.009191616 0.4404099 -0.02015771 -0.05450731
## 2009-02-02 0.02839804 0.073164614 0.4960102 0.05844322 0.02436688
## d_22 d_23 d_25 d_26 d_28
## 2008-12-29 -0.32997381 -0.04321831 -0.308494915 -0.27304142 -0.319469993
## 2009-01-05 -0.24930346 0.01694022 -0.237538275 -0.19927446 -0.248989923
## 2009-01-12 -0.09323119 0.13976207 -0.091931001 -0.04998092 -0.103850614
## 2009-01-19 -0.11675043 0.11389645 -0.123396553 -0.07999985 -0.135764867
## 2009-01-26 -0.05583703 0.15848662 -0.070898689 -0.02515271 -0.083685784
## 2009-02-02 0.02675053 0.22160476 0.003724671 0.05193142 -0.009442349
## d_29 d_3 d_31 d_32 d_33
## 2008-12-29 -0.5856532 -0.246254713 -0.04435051 -0.39912964 -0.42176903
## 2009-01-05 -0.4916583 -0.172599084 0.00969803 -0.32176695 -0.34102991
## 2009-01-12 -0.3237150 -0.036523860 0.13994159 -0.16891432 -0.18209195
## 2009-01-19 -0.3331191 -0.049681662 0.09095306 -0.19558441 -0.20916787
## 2009-01-26 -0.2603395 0.006768922 0.12843233 -0.13758125 -0.14855719
## 2009-02-02 -0.1674179 0.080647053 0.18987106 -0.05762385 -0.06524941
## d_34 d_35 d_36 d_37 d_38
## 2008-12-29 -0.28112122 -0.49119746 -0.212012955 -0.224826755 -0.31638715
## 2009-01-05 -0.20327378 -0.40254675 -0.145301379 -0.158012155 -0.23777916
## 2009-01-12 -0.04717529 -0.24312537 -0.002907525 -0.015967457 -0.08370372
## 2009-01-19 -0.07697380 -0.25332356 -0.039605865 -0.051985985 -0.10920133
## 2009-01-26 -0.01890415 -0.18481517 0.009037761 -0.003193387 -0.05010388
## 2009-02-02 0.06209836 -0.09693262 0.080504014 0.068247914 0.03084629
## d_39 d_4 d_40 d_41 d_42
## 2008-12-29 -0.29213926 -0.29973410 -0.53222152 -0.32809583 -0.378155279
## 2009-01-05 -0.22279843 -0.21398232 -0.43357580 -0.25764133 -0.291276109
## 2009-01-12 -0.07753872 -0.02041942 -0.22735179 -0.11202246 -0.096606178
## 2009-01-19 -0.11212922 -0.08121781 -0.27601026 -0.14461394 -0.156343111
## 2009-01-26 -0.06121204 -0.02005713 -0.20351918 -0.09262282 -0.094191751
## 2009-02-02 0.01241783 0.07430993 -0.09887297 -0.01827976 0.001074077
## d_43 d_44 d_45 d_46 d_47
## 2008-12-29 -0.34548522 -0.360236872 -0.3829384 -0.27244354 -0.27882706
## 2009-01-05 -0.27163435 -0.280152836 -0.3073032 -0.19930607 -0.20437111
## 2009-01-12 -0.12223012 -0.124656291 -0.1565973 -0.05155735 -0.05160286
## 2009-01-19 -0.15220666 -0.148727555 -0.1843111 -0.08096961 -0.08459452
## 2009-01-26 -0.09728946 -0.088329369 -0.1277676 -0.02655536 -0.02950509
## 2009-02-02 -0.02013171 -0.006209217 -0.0492945 0.04980657 0.04879370
## d_48 d_49 d_5 d_50 d_51
## 2008-12-29 -0.27181160 -0.26771743 -0.374734061 -0.319956942 -0.45751696
## 2009-01-05 -0.19942577 -0.19857965 -0.290554370 -0.246515366 -0.37468260
## 2009-01-12 -0.05241511 -0.06694070 -0.130723757 -0.098468031 -0.22097254
## 2009-01-19 -0.08253505 -0.08435216 -0.151343153 -0.127593964 -0.23664695
## 2009-01-26 -0.02878130 -0.03187159 -0.087386499 -0.072912494 -0.17324960
## 2009-02-02 0.04698142 0.03840488 -0.001926981 0.003691875 -0.09000388
## d_52 d_53 d_54 d_55 d_56
## 2008-12-29 -0.27572227 -0.19067775 -0.30915499 -0.5262448 -0.35429451
## 2009-01-05 -0.20163791 -0.11506116 -0.23279978 -0.4331306 -0.27198606
## 2009-01-12 -0.05245473 0.06854956 -0.07816656 -0.2323382 -0.11399290
## 2009-01-19 -0.08162863 -0.00179134 -0.10937003 -0.2862047 -0.13637411
## 2009-01-26 -0.02644782 0.05046316 -0.05261164 -0.2185743 -0.07406179
## 2009-02-02 0.05078899 0.13675037 0.02720125 -0.1183379 0.00990596
## d_57 d_58 d_6 d_60 d_61
## 2008-12-29 -0.5090624 -0.32387430 -0.31846047 -0.32444892 -0.27866853
## 2009-01-05 -0.4211834 -0.24304196 -0.24506704 -0.24709923 -0.20002604
## 2009-01-12 -0.2592455 -0.08802352 -0.09614028 -0.09148950 -0.04715790
## 2009-01-19 -0.2744079 -0.10982048 -0.12651089 -0.12175665 -0.07101665
## 2009-01-26 -0.2070025 -0.04860732 -0.07199198 -0.06412439 -0.01172779
## 2009-02-02 -0.1189566 0.03382096 0.00479437 0.01648131 0.06895473
## d_62 d_63 d_64 d_65 d_66
## 2008-12-29 -0.35613912 -0.37540400 -0.43521285 -0.43792510 -0.27135513
## 2009-01-05 -0.27299408 -0.29561931 -0.35699482 -0.34845599 -0.20017004
## 2009-01-12 -0.11449180 -0.14490384 -0.20375260 -0.18371556 -0.05309936
## 2009-01-19 -0.13568106 -0.16344958 -0.22903455 -0.19898720 -0.08595345
## 2009-01-26 -0.07259305 -0.10273209 -0.17022132 -0.13034567 -0.03341564
## 2009-02-02 0.01196736 -0.02191759 -0.08968912 -0.04073695 0.04168449
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59
## 2008-12-29 -0.19562538 -0.370729639 -0.26258597 -0.28860483 NA NA NA NA
## 2009-01-05 -0.12208967 -0.286260771 -0.18629692 -0.21223366 NA NA NA NA
## 2009-01-12 0.05947775 -0.130945720 -0.03573972 -0.06487009 NA NA NA NA
## 2009-01-19 -0.01282235 -0.145081208 -0.06181428 -0.08662973 NA NA NA NA
## 2009-01-26 0.03760360 -0.080247562 -0.00459347 -0.02891185 NA NA NA NA
## 2009-02-02 0.12223192 0.004301193 0.07421288 0.04918136 NA NA NA NA
## d_9
## 2008-12-29 NA
## 2009-01-05 NA
## 2009-01-12 NA
## 2009-01-19 NA
## 2009-01-26 NA
## 2009-02-02 NAtail(pred.A.cv$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 1.0150622 1.0382420 1.0296746 1.0688803 0.9965445 0.9528974
## 2020-11-09 NA 0.6131506 0.6328900 0.6109485 0.6540798 0.5659699 0.5474843
## 2020-11-16 NA 0.9610034 1.0004710 0.9729153 1.0199199 0.9413175 0.8996281
## 2020-11-23 NA 0.8772879 0.9082165 0.8907819 0.9269566 0.8506757 0.7985353
## 2020-11-30 NA 1.0235842 1.0585015 1.0497072 1.0851006 1.0173697 0.9517047
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2
## 2020-11-02 1.0804114 1.1536438 1.1680318 1.0220347 1.0986978 1.3051833
## 2020-11-09 0.6764669 0.7304555 0.7529365 0.6150522 0.6973173 0.8564288
## 2020-11-16 1.0460786 1.1273576 1.1093427 0.9744087 1.0664945 1.2412267
## 2020-11-23 0.9450127 1.0175735 1.0160217 0.8856271 0.9711145 1.1597817
## 2020-11-30 1.0858757 1.2088914 1.1650578 1.0363943 1.1146498 1.3210208
## 2020-12-07 NA NA NA NA NA NA
## d_20 d_21 d_22 d_23 d_25 d_26
## 2020-11-02 1.0285957 1.0365695 1.0708317 1.0895563 0.9535186 1.0454334
## 2020-11-09 0.6213195 0.6251947 0.6622079 0.6833695 0.5476990 0.6133030
## 2020-11-16 0.9962454 0.9876766 1.0168746 1.0645528 0.9155780 0.9862652
## 2020-11-23 0.9052317 0.8990545 0.9306350 0.9657703 0.8207370 0.9015860
## 2020-11-30 1.0750396 1.0501848 1.0797193 1.1212032 0.9763268 1.0626299
## 2020-12-07 NA NA NA NA NA NA
## d_28 d_29 d_3 d_31 d_32 d_33
## 2020-11-02 0.9161651 1.0337061 1.1157395 0.9396360 0.9541425 0.9523659
## 2020-11-09 0.5196546 0.6279831 0.6933270 0.5406458 0.5485554 0.5574240
## 2020-11-16 0.8770143 0.9779460 1.0643686 0.9160772 0.9081766 0.9025589
## 2020-11-23 0.7859653 0.8833830 0.9695286 0.8271749 0.8169419 0.8182352
## 2020-11-30 0.9375273 1.0281945 1.1289236 0.9803210 0.9674177 0.9614071
## 2020-12-07 NA NA NA NA NA NA
## d_34 d_35 d_36 d_37 d_38 d_39
## 2020-11-02 1.0803700 1.0705840 0.9696009 0.9661821 1.0635672 0.9032826
## 2020-11-09 0.6680280 0.6534341 0.5664519 0.5601490 0.6515349 0.5181359
## 2020-11-16 1.0341767 1.0155012 0.9294754 0.9279890 1.0164734 0.8636166
## 2020-11-23 0.9427047 0.9178723 0.8418020 0.8379460 0.9210879 0.7803393
## 2020-11-30 1.0982584 1.0725128 0.9945305 0.9918441 1.0769120 0.9244404
## 2020-12-07 NA NA NA NA NA NA
## d_4 d_40 d_41 d_42 d_43 d_44
## 2020-11-02 1.0434203 0.8566342 0.9275682 0.9986371 0.9570589 1.0351823
## 2020-11-09 0.6315746 0.4780776 0.5178459 0.5906635 0.5508841 0.6235685
## 2020-11-16 0.9848015 0.7916788 0.8884220 0.9476964 0.9154993 0.9808281
## 2020-11-23 0.9139090 0.6923339 0.7918549 0.8708545 0.8210061 0.8943801
## 2020-11-30 1.0573747 0.8180439 0.9512450 1.0117757 0.9757378 1.0448760
## 2020-12-07 NA NA NA NA NA NA
## d_45 d_46 d_47 d_48 d_49 d_5
## 2020-11-02 0.9437724 0.9964143 1.0125199 0.9969173 1.0191359 1.0208804
## 2020-11-09 0.5391465 0.5978956 0.6169192 0.5959013 0.6051918 0.6261924
## 2020-11-16 0.9006534 0.9473731 0.9830464 0.9521620 0.9803016 0.9608471
## 2020-11-23 0.8074426 0.8651166 0.8934628 0.8659315 0.8832129 0.8734667
## 2020-11-30 0.9597959 1.0129488 1.0285131 1.0150767 1.0375160 1.0161722
## 2020-12-07 NA NA NA NA NA NA
## d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 0.9448261 0.9990481 1.0295018 1.1381029 1.0084533 0.9396546
## 2020-11-09 0.5560352 0.5963697 0.6156387 0.6914741 0.6008224 0.5421431
## 2020-11-16 0.8995586 0.9609164 0.9785194 1.0678106 0.9713484 0.8915966
## 2020-11-23 0.8168895 0.8646994 0.8921744 0.9928311 0.8825779 0.8102778
## 2020-11-30 0.9606547 1.0048405 1.0459417 1.1407781 1.0457167 0.9464136
## 2020-12-07 NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61
## 2020-11-02 1.0344882 1.0033416 1.1211895 0.9775367 0.9970461 1.0676910
## 2020-11-09 0.6331909 0.6044010 0.7026047 0.5683030 0.5912209 0.6601302
## 2020-11-16 0.9741131 0.9514381 1.0680028 0.9333517 0.9465716 1.0080392
## 2020-11-23 0.8913156 0.8599500 0.9774114 0.8433598 0.8560768 0.9249415
## 2020-11-30 1.0374271 1.0040557 1.1261933 0.9954516 1.0069950 1.0754843
## 2020-12-07 NA NA NA NA NA NA
## d_62 d_63 d_64 d_65 d_66 d_67
## 2020-11-02 1.0473612 1.0332122 0.9144571 1.1152067 0.9806742 0.9883702
## 2020-11-09 0.6458638 0.6242988 0.5192354 0.6730892 0.5727237 0.5843624
## 2020-11-16 0.9863050 0.9817355 0.8683077 1.0385928 0.9371660 0.9611449
## 2020-11-23 0.9035327 0.8913614 0.7798741 0.9551312 0.8456867 0.8777696
## 2020-11-30 1.0494718 1.0418505 0.9257006 1.1170440 1.0015214 1.0163362
## 2020-12-07 NA NA NA NA NA NA
## d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.1616029 1.0505711 1.0690504 NA NA NA NA NA
## 2020-11-09 0.7352913 0.6397712 0.6581088 NA NA NA NA NA
## 2020-11-16 1.0998327 0.9929101 1.0187577 NA NA NA NA NA
## 2020-11-23 1.0077940 0.9066352 0.9293125 NA NA NA NA NA
## 2020-11-30 1.1584806 1.0596425 1.0815680 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA NA NAhead(pred.A.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 0.7072519 0.8107323 0.8580920 0.8138605 0.7877827 0.6602204
## 2009-01-05 NA 0.7664933 0.8708468 0.9193588 0.8731730 0.8421988 0.7143733
## 2009-01-12 NA 0.8960555 1.0076308 1.0619536 0.9972429 0.9585195 0.8227071
## 2009-01-19 NA 0.8747138 0.9773225 1.0265238 0.9812117 0.9399631 0.8159760
## 2009-01-26 NA 0.9294706 1.0305796 1.0799825 1.0352766 0.9886653 0.8673813
## 2009-02-02 NA 1.0094390 1.1108872 1.1620660 1.1118419 1.0587820 0.9378538
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2008-12-29 0.6754515 1.280104 0.6775480 0.7552597 0.7767860 1.381085 0.7792113
## 2009-01-05 0.7361991 1.338311 0.7414268 0.8145943 0.8409609 1.446645 0.8402981
## 2009-01-12 0.8643781 1.490123 0.8708222 0.9475350 0.9845313 1.635079 0.9768541
## 2009-01-19 0.8496685 1.430878 0.8631942 0.9211346 0.9570371 1.546461 0.9500682
## 2009-01-26 0.9075086 1.475776 0.9256275 0.9747172 1.0156684 1.594985 1.0052790
## 2009-02-02 0.9897528 1.552557 1.0118845 1.0544713 1.1027668 1.685815 1.0873124
## d_21 d_22 d_23 d_25 d_26 d_28
## 2008-12-29 0.7498693 0.7392707 0.9841391 0.7533963 0.7825803 0.7451729
## 2009-01-05 0.8093134 0.8010269 1.0447848 0.8086894 0.8421208 0.7994814
## 2009-01-12 0.9408188 0.9359662 1.1809465 0.9353352 0.9773351 0.9242526
## 2009-01-19 0.9166182 0.9139083 1.1504822 0.9062722 0.9481198 0.8951321
## 2009-01-26 0.9706181 0.9710448 1.2026769 0.9550422 1.0013025 0.9429077
## 2009-02-02 1.0502304 1.0544187 1.2808092 1.0289702 1.0813070 1.0155105
## d_29 d_3 d_31 d_32 d_33 d_34
## 2008-12-29 0.5708173 0.8033021 0.9815299 0.6882465 0.6726741 0.7742608
## 2009-01-05 0.6270192 0.8643925 1.0358470 0.7434297 0.7291175 0.8368050
## 2009-01-12 0.7416252 0.9900799 1.1797459 0.8660299 0.8545948 0.9780337
## 2009-01-19 0.7346358 0.9768746 1.1231879 0.8430892 0.8316633 0.9492021
## 2009-01-26 0.7900537 1.0333770 1.1659492 0.8933074 0.8835401 1.0058514
## 2009-02-02 0.8669489 1.1124161 1.2397170 0.9675556 0.9602211 1.0906306
## d_35 d_36 d_37 d_38 d_39 d_4
## 2008-12-29 0.6277744 0.8298656 0.8200733 0.7476160 0.7661084 0.7608680
## 2009-01-05 0.6858281 0.8869066 0.8764566 0.8085657 0.8209632 0.8287259
## 2009-01-12 0.8042165 1.0224185 1.0099430 0.9430601 0.9491569 1.0054282
## 2009-01-19 0.7959346 0.9854036 0.9739792 0.9191559 0.9167584 0.9458925
## 2009-01-26 0.8522667 1.0343720 1.0224784 0.9749714 0.9645353 1.0053501
## 2009-02-02 0.9304628 1.1108722 1.0980256 1.0570561 1.0381392 1.1046685
## d_40 d_41 d_42 d_43 d_44 d_45
## 2008-12-29 0.6030333 0.7396112 0.7034795 0.7261753 0.7172333 0.7001413
## 2009-01-05 0.6653386 0.7933449 0.7670835 0.7816498 0.7766930 0.7549082
## 2009-01-12 0.8174881 0.9174464 0.9316731 0.9074183 0.9070107 0.8774495
## 2009-01-19 0.7784747 0.8878138 0.8774356 0.8804653 0.8851465 0.8532605
## 2009-01-26 0.8368368 0.9350082 0.9335145 0.9300353 0.9400003 0.9027183
## 2009-02-02 0.9290083 1.0070104 1.0266586 1.0045199 1.0202318 0.9762576
## d_46 d_47 d_48 d_49 d_5 d_50
## 2008-12-29 0.7810528 0.7760392 0.7815465 0.7862448 0.7053745 0.7448101
## 2009-01-05 0.8402060 0.8358872 0.8401054 0.8422243 0.7671807 0.8014622
## 2009-01-12 0.9738708 0.9737130 0.9730358 0.9604182 0.8999948 0.9292408
## 2009-01-19 0.9455498 0.9419960 0.9440708 0.9435863 0.8815045 0.9024761
## 2009-01-26 0.9983449 0.9952448 0.9961251 0.9942084 0.9396167 0.9531209
## 2009-02-02 1.0774966 1.0762163 1.0744568 1.0664039 1.0233536 1.0289364
## d_51 d_52 d_53 d_54 d_55 d_56
## 2008-12-29 0.6492782 0.7793796 0.8485392 0.7528568 0.6066483 0.7199404
## 2009-01-05 0.7052069 0.8390426 0.9148962 0.8124587 0.6656348 0.7815590
## 2009-01-12 0.8222310 0.9737571 1.0989800 0.9481881 0.8134220 0.9151792
## 2009-01-19 0.8093194 0.9455316 1.0240856 0.9189443 0.7705789 0.8947990
## 2009-01-26 0.8621809 0.9989755 1.0788073 0.9725118 0.8243324 0.9522206
## 2009-02-02 0.9369321 1.0790216 1.1758435 1.0532272 0.9111001 1.0355349
## d_57 d_58 d_6 d_60 d_61 d_62
## 2008-12-29 0.6162545 0.7416282 0.7478314 0.7414303 0.7759234 0.7201783
## 2009-01-05 0.6728027 0.8039954 0.8044275 0.8009236 0.8393346 0.7822731
## 2009-01-12 0.7910123 0.9387371 0.9332473 0.9356392 0.9778938 0.9162770
## 2009-01-19 0.7790584 0.9184371 0.9050316 0.9076319 0.9547765 0.8967702
## 2009-01-26 0.8333368 0.9763617 0.9554836 0.9613797 1.0130416 0.9549095
## 2009-02-02 0.9099970 1.0602064 1.0315200 1.0419970 1.0981175 1.0389457
## d_63 d_64 d_65 d_66 d_67 d_68
## 2008-12-29 0.7048424 0.6644819 0.6620570 0.7821979 0.8443513 0.7081448
## 2009-01-05 0.7632263 0.7183124 0.7238507 0.8397521 0.9084884 0.7704025
## 2009-01-12 0.8872174 0.8370335 0.8533059 0.9726394 1.0890553 0.8996881
## 2009-01-19 0.8707815 0.8159405 0.8402252 0.9410721 1.0128510 0.8869241
## 2009-01-26 0.9251747 0.8651966 0.8997944 0.9917220 1.0650232 0.9462124
## 2009-02-02 1.0029461 0.9376080 0.9840334 1.0689711 1.1588954 1.0295899
## d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.7885031 0.7687558 NA NA NA NA NA
## 2009-01-05 0.8509374 0.8295971 NA NA NA NA NA
## 2009-01-12 0.9891236 0.9611435 NA NA NA NA NA
## 2009-01-19 0.9636032 0.9403112 NA NA NA NA NA
## 2009-01-26 1.0202949 0.9960554 NA NA NA NA NA
## 2009-02-02 1.1039068 1.0768507 NA NA NA NA NAtail(pred.A.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 2.781599 2.838884 2.819998 2.931073 2.724446 2.610302
## 2020-11-09 NA 1.860912 1.892757 1.855087 1.935722 1.771103 1.740139
## 2020-11-16 NA 2.634997 2.733535 2.664019 2.790584 2.577655 2.474502
## 2020-11-23 NA 2.423315 2.492598 2.453837 2.542715 2.354157 2.236453
## 2020-11-30 NA 2.805015 2.896775 2.876397 2.978239 2.781061 2.606521
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2020-11-02 2.958878 3.183435 3.235531 2.795816 3.014295 3.709308 2.811627
## 2020-11-09 1.975516 2.084823 2.136237 1.860942 2.017667 2.367855 1.870878
## 2020-11-16 2.858830 3.100352 3.050820 2.665503 2.918542 3.478839 2.721738
## 2020-11-23 2.583961 2.777846 2.778904 2.438966 2.652962 3.206522 2.484836
## 2020-11-30 2.974775 3.363389 3.225430 2.835765 3.062372 3.767352 2.944605
## 2020-12-07 NA NA NA NA NA NA NA
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2020-11-02 2.835976 2.942732 2.986667 2.609418 2.864462 2.516620 2.829909
## 2020-11-09 1.879464 1.955429 1.989499 1.738933 1.859197 1.692778 1.886074
## 2020-11-16 2.700533 2.787634 2.912454 2.512108 2.699383 2.419842 2.676315
## 2020-11-23 2.471466 2.557135 2.638360 2.284754 2.480047 2.209201 2.434795
## 2020-11-30 2.874646 2.968091 3.081899 2.669337 2.913247 2.570692 2.814154
## 2020-12-07 NA NA NA NA NA NA NA
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2020-11-02 3.069518 2.572590 2.611916 2.612236 2.961643 2.932664 2.652622
## 2020-11-09 2.011780 1.726117 1.740963 1.759838 1.960812 1.932296 1.772411
## 2020-11-16 2.915356 2.512473 2.494312 2.485125 2.827726 2.775232 2.548023
## 2020-11-23 2.651423 2.298679 2.276733 2.284101 2.580475 2.517011 2.334056
## 2020-11-30 3.109453 2.679035 2.646367 2.635632 3.014715 2.937874 2.719115
## 2020-12-07 NA NA NA NA NA NA NA
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2020-11-02 2.644205 2.911453 2.488470 2.863516 2.372914 2.542926 2.736027
## 2020-11-09 1.761663 1.928150 1.692952 1.896721 1.624961 1.687948 1.819306
## 2020-11-16 2.544758 2.777227 2.391491 2.700107 2.223367 2.444956 2.599777
## 2020-11-23 2.325518 2.524468 2.200336 2.515193 2.013004 2.219784 2.407369
## 2020-11-30 2.712311 2.950057 2.541329 2.903091 2.282562 2.603249 2.771575
## 2020-12-07 NA NA NA NA NA NA NA
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2020-11-02 2.618288 2.838952 2.586647 2.729897 2.766253 2.728895 2.784979
## 2020-11-09 1.744184 1.880836 1.725743 1.832553 1.862369 1.827311 1.840814
## 2020-11-16 2.511438 2.688256 2.477133 2.599086 2.685701 2.609290 2.678475
## 2020-11-23 2.284907 2.465460 2.256560 2.393797 2.455501 2.393664 2.430515
## 2020-11-30 2.667195 2.865725 2.627821 2.775124 2.810490 2.778615 2.835907
## 2020-12-07 NA NA NA NA NA NA NA
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 2.798606 2.593982 2.729173 2.818667 3.138277 2.755226 2.580497
## 2020-11-09 1.885864 1.758336 1.824430 1.863250 2.007652 1.832763 1.733933
## 2020-11-16 2.635317 2.479021 2.626816 2.678187 2.924839 2.654659 2.459072
## 2020-11-23 2.414743 2.282275 2.385770 2.456520 2.713420 2.429094 2.266906
## 2020-11-30 2.785075 2.635097 2.744604 2.864727 3.145953 2.859466 2.597403
## 2020-12-07 NA NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 2.837368 2.746733 3.084215 2.677513 2.726781 2.930210 2.878161
## 2020-11-09 1.899389 1.843097 2.029294 1.778107 1.817121 1.949337 1.926201
## 2020-11-16 2.670904 2.607692 2.924332 2.561299 2.592363 2.760409 2.707182
## 2020-11-23 2.458594 2.379670 2.671016 2.340715 2.368005 2.540259 2.491962
## 2020-11-30 2.845325 2.748500 3.099462 2.725073 2.753696 2.952927 2.883361
## 2020-12-07 NA NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2020-11-02 2.828102 2.515289 3.068189 2.681911 2.701232 3.211837 2.878055
## 2020-11-09 1.878813 1.693992 1.971730 1.783417 1.803304 2.096940 1.908448
## 2020-11-16 2.685954 2.401511 2.841604 2.567506 2.628308 3.019160 2.716679
## 2020-11-23 2.453778 2.198148 2.613969 2.342984 2.417937 2.753594 2.492085
## 2020-11-30 2.852203 2.543146 3.073304 2.738022 2.777195 3.201332 2.904075
## 2020-12-07 NA NA NA NA NA NA NA
## d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.930663 NA NA NA NA NA
## 2020-11-09 1.943003 NA NA NA NA NA
## 2020-11-16 2.786657 NA NA NA NA NA
## 2020-11-23 2.548142 NA NA NA NA NA
## 2020-11-30 2.967125 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NAsummary(pred.A.cv)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX
## obs 0.16890294 0.16890294 0.10978414
## average 0.09626864 0.09626864 0.05356057
##
## R2:
## EX.mu EX.mu.beta EX
## obs 0.5764898 0.5764898 0.8210760
## average 0.5882971 0.5882971 0.8725603
##
## Coverage of 95% prediction intervals:
## EX
## obs 0.9405063
## average 0.8412698summary(pred.A.cv.log)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.2193848 0.2193848 0.15105214 0.15107290
## average 0.1077456 0.1077456 0.06196514 0.06244725
##
## R2:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.6069412 0.6069412 0.8136634 0.8136122
## average 0.6550145 0.6550145 0.8858972 0.8841148
##
## Coverage of 95% prediction intervals:
## EX.pred
## obs 0.9405063
## average 0.8888889par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModA.jpeg"),
width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()## quartz_off_screen
## 23.1.5 Predictions for 2009-2020
What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites
x.A <- coef(est.denver.model.A, pars = "cov")$par
x.A## [1] -7.789155 -7.377014 12.446866 -3.007835 -4.861449E.A <- predict(denver.model.A, est.denver.model.A, LTA = T,
transform="unbiased", pred.var=FALSE)
print(E.A)## Prediction for STmodel.
##
## Regression parameters:
## 0 Spatio-temporal covariate(s).
## 26 beta-fields regression parameters in x$pars.
##
## Regression parameters are assumed to be known and
## prediction variances do NOT include
## uncertainties in regression parameters.
##
## Prediction of beta-fields, (x$beta):
## List of 2
## $ mu: num [1:69, 1:2] 0.2513 0.0254 0.0394 0.0729 0.0457 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX: num [1:69, 1:2] 0.2601 0.0396 0.0359 0.063 0.0316 ...
## ..- attr(*, "dimnames")=List of 2
##
## Predictions for log-Gaussian field of type: unbiased
##
## Predictions for 622 times at 69 locations.
## List of 5
## $ EX.mu : num [1:622, 1:69] 1.2 1.26 1.43 1.36 1.4 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.mu.beta: num [1:622, 1:69] 1.26 1.32 1.48 1.4 1.45 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX : num [1:622, 1:69] 1.29 1.35 1.52 1.44 1.48 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.pred : num [1:622, 1:69] 1.29 1.35 1.53 1.45 1.49 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.EX : num [1:622, 1:69] 0.23 0.275 0.395 0.34 0.37 ...
## ..- attr(*, "dimnames")=List of 2
##
## Variances have been computed.
## List of 2
## $ log.VX : num [1:622, 1:69] 0.0499 0.0498 0.0498 0.0497 0.0497 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.VX.pred: num [1:622, 1:69] 0.0577 0.0576 0.0575 0.0574 0.0574 ...
## ..- attr(*, "dimnames")=List of 2
##
## Mean squared prediciton errors NOT been computed.
##
## 69 temporal averages have been compute.
## List of 1
## $ LTA:'data.frame': 69 obs. of 4 variables:week_preds <- as.data.frame(E.A$EX) %>%
mutate(week = as.Date(rownames(E.A$EX)),
year = year(as.Date(rownames(E.A$EX))))
week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]
box_preds <- week_preds %>%
select(week, all_of(week_sites)) %>%
#filter(week %in% week_weeks) %>%
mutate(month = month(week),
year = year(week)) %>%
filter(year %in% c(2009:2020))
box_data <- box_preds %>%
pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)## week month year location
## Min. :2009-01-05 Min. : 1.000 Min. :2009 Length:42228
## 1st Qu.:2012-01-09 1st Qu.: 4.000 1st Qu.:2012 Class :character
## Median :2014-12-29 Median : 7.000 Median :2014 Mode :character
## Mean :2014-12-28 Mean : 6.494 Mean :2014
## 3rd Qu.:2017-12-18 3rd Qu.: 9.000 3rd Qu.:2017
## Max. :2020-12-07 Max. :12.000 Max. :2020
##
## pred
## Min. :0.3178
## 1st Qu.:1.1234
## Median :1.3756
## Mean :1.4203
## 3rd Qu.:1.6490
## Max. :3.9057
## NA's :861#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
group_by(month) %>%
summarize(median = median(pred, na.rm=T)) %>%
arrange(desc(median))
box_summarypred_box_plot <- ggplot(box_data) +
geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
show.legend = F) +
#scale_color_viridis(name = "Month", discrete = T) +
facet_wrap(. ~ year, scales = "free") +
xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
# theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot
3.2 Model B: 2 temporal trend (basis) functions
3.2.1 Create the model object
For this version of the model, use iid for cov.beta (beta, beta1, beta2) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.
denver.data.B <- denver.data2.2
names(denver.data.B$covars)## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "high_int_100"
## [9] "pop_ct_1000" "dist_m_cafos" "dist_m_og_wells" "dist_m_wwtp"
## [13] "aadt_100" "aadt_250" "aadt_1000" "x"
## [17] "y" "type"LUR.B <- list(covar_fun, covar_fun, covar_fun)
cov.beta.B <- list(covf = c("iid", "iid", "iid"), nugget = T)
cov.nu.B <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.B <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")
denver.model.B <- createSTmodel(denver.data.B, LUR = LUR.B,
ST = c("bc_st_no2", "bc_st_smk"),
cov.beta = cov.beta.B, cov.nu = cov.nu.B,
locations = locations.B)## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.denver.model.B## STmodel-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## Models for the beta-fields are:
## $const
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## high_int_100 + pop_ct_1000 + dist_m_cafos + dist_m_og_wells +
## dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## $V1
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## high_int_100 + pop_ct_1000 + dist_m_cafos + dist_m_og_wells +
## dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## $V2
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## high_int_100 + pop_ct_1000 + dist_m_cafos + dist_m_og_wells +
## dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## Covariance model for the beta-field(s):
## Covariance type(s): iid, iid, iid
## Nugget: Yes, Yes, Yes
## Covariance model for the nu-field(s):
## Covariance type: exp
## Nugget: ~1
## Random effect: No
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-063.2.2 Estimate model parameters
- nugget values can range from -1 to -6
- range values can range from 0 to 4
- try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
- make sure the starting values are pretty different
names <- loglikeSTnames(denver.model.B, all=FALSE)
names## [1] "log.nugget.const.iid" "log.nugget.V1.iid"
## [3] "log.nugget.V2.iid" "nu.log.range.exp"
## [5] "nu.log.sill.exp" "nu.log.nugget.(Intercept).exp"# x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
# c(-1, -1, -1, 0, -1, -1),
# c(-1, -1, -1, 0, -5, -1),
# c(-5, -5, -5, 0, -1, -5),
# c(-5, -5, -5, 0, -5, -5),
# c(-1, -1, -1, 2, -1, -1),
# c(-1, -1, -1, 2, -5, -1),
# c(-5, -5, -5, 2, -1, -5),
# c(-5, -5, -5, 2, -5, -5),
# c(-1, -1, -1, 4, -1, -1),
# c(-1, -1, -1, 4, -5, -1),
# c(-5, -5, -5, 4, -1, -5),
# c(-5, -5, -5, 4, -5, -5))
x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
c(-10, -5, -5, 4, -3, -5),
c(-10, -5, -5, 4, -5, -5),
c(-10, -5, -5, 6, -3, -5),
c(-10, -5, -5, 6, -5, -5),
c(-10, -5, -5, 8, -3, -5),
c(-10, -5, -5, 8, -5, -5),
c(-10, -5, -5, 10, -3, -5),
c(-10, -5, -5, 10, -5, -5),
c(-12, -5, -5, 8, -5, -5),
c(-12, -5, -5, 10, -5, -5),
c(-12, -5, -5, 12, -5, -5))
rownames(x.init.B) <- loglikeSTnames(denver.model.B, all=FALSE)
x.init.B## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## log.nugget.const.iid 0 -10 -10 -10 -10 -10 -10 -10 -10
## log.nugget.V1.iid 0 -5 -5 -5 -5 -5 -5 -5 -5
## log.nugget.V2.iid 0 -5 -5 -5 -5 -5 -5 -5 -5
## nu.log.range.exp 0 4 4 6 6 8 8 10 10
## nu.log.sill.exp 0 -3 -5 -3 -5 -3 -5 -3 -5
## nu.log.nugget.(Intercept).exp 0 -5 -5 -5 -5 -5 -5 -5 -5
## [,10] [,11] [,12]
## log.nugget.const.iid -12 -12 -12
## log.nugget.V1.iid -5 -5 -5
## log.nugget.V2.iid -5 -5 -5
## nu.log.range.exp 8 10 12
## nu.log.sill.exp -5 -5 -5
## nu.log.nugget.(Intercept).exp -5 -5 -5#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.B <- estimate.STmodel(denver.model.B, x.init.B)## Optimisation using starting value 1/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= 507.58 |proj g|= 15
## At iterate 10 f = -1392.9 |proj g|= 2.3574
## At iterate 20 f = -1393.2 |proj g|= 0.21723
##
## iterations 26
## function evaluations 33
## segments explored during Cauchy searches 31
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 2
## norm of the final projected gradient 0.0938135
## final function value -1393.21
##
## F = -1393.21
## final value -1393.211978
## converged## Optimisation using starting value 2/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1197.6 |proj g|= 12
## At iterate 10 f = -1392.8 |proj g|= 1.258
## At iterate 20 f = -1394.1 |proj g|= 11.199
## ys=-3.157e+00 -gs= 2.729e+00, BFGS update SKIPPED
## At iterate 30 f = -1619.2 |proj g|= 10.221
## At iterate 40 f = -1620.7 |proj g|= 0.53425
## At iterate 50 f = -1632.5 |proj g|= 1.2872
## At iterate 60 f = -1635.2 |proj g|= 2.5951
##
## iterations 65
## function evaluations 111
## segments explored during Cauchy searches 68
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00710589
## final function value -1635.44
##
## F = -1635.44
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1635.442472
## converged## Optimisation using starting value 3/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1237.3 |proj g|= 20
## At iterate 10 f = -1391.3 |proj g|= 2.4494
## At iterate 20 f = -1393.1 |proj g|= 0.64933
## ys=-2.727e+01 -gs= 8.698e-01, BFGS update SKIPPED
## At iterate 30 f = -1615.6 |proj g|= 7.7933
## At iterate 40 f = -1620.7 |proj g|= 0.024487
## At iterate 50 f = -1632.1 |proj g|= 1.5927
## ys=-1.790e-01 -gs= 2.422e-02, BFGS update SKIPPED
## At iterate 60 f = -1632.3 |proj g|= 0.38934
## At iterate 70 f = -1634.9 |proj g|= 3.8621
## At iterate 80 f = -1635.4 |proj g|= 0.031101
##
## iterations 82
## function evaluations 131
## segments explored during Cauchy searches 86
## BFGS updates skipped 2
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0196387
## final function value -1635.45
##
## F = -1635.45
## final value -1635.445721
## converged## Optimisation using starting value 4/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1202.7 |proj g|= 12
## At iterate 10 f = -1624.9 |proj g|= 10.136
## At iterate 20 f = -1635.1 |proj g|= 1.4066
## At iterate 30 f = -1635.4 |proj g|= 0.67206
## At iterate 40 f = -1635.4 |proj g|= 0.46818
## At iterate 50 f = -1635.4 |proj g|= 0.80816
## At iterate 60 f = -1635.4 |proj g|= 0.036044
## At iterate 70 f = -1635.4 |proj g|= 0.021507
##
## iterations 71
## function evaluations 89
## segments explored during Cauchy searches 75
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0208085
## final function value -1635.45
##
## F = -1635.45
## final value -1635.445728
## converged## Optimisation using starting value 5/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1244.8 |proj g|= 20
## At iterate 10 f = -1612.9 |proj g|= 3.223
## At iterate 20 f = -1623.7 |proj g|= 10.307
## At iterate 30 f = -1625.6 |proj g|= 8.858
## At iterate 40 f = -1628.2 |proj g|= 10.126
##
## iterations 48
## function evaluations 68
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0262019
## final function value -1629.64
##
## F = -1629.64
## final value -1629.644361
## converged## Optimisation using starting value 6/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1262.8 |proj g|= 12
## At iterate 10 f = -1635.2 |proj g|= 2.8922
## At iterate 20 f = -1635.4 |proj g|= 0.14032
## At iterate 30 f = -1635.4 |proj g|= 0.99514
## At iterate 40 f = -1635.4 |proj g|= 0.039377
## At iterate 50 f = -1635.4 |proj g|= 0.039755
##
## iterations 50
## function evaluations 59
## segments explored during Cauchy searches 55
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0397554
## final function value -1635.45
##
## F = -1635.45
## final value -1635.445653
## converged## Optimisation using starting value 7/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1297.4 |proj g|= 20
## At iterate 10 f = -1634.6 |proj g|= 2.5526
## At iterate 20 f = -1635.4 |proj g|= 0.31687
## At iterate 30 f = -1635.4 |proj g|= 0.27392
## At iterate 40 f = -1635.4 |proj g|= 0.05911
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 45
## function evaluations 77
## segments explored during Cauchy searches 51
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0267835
## final function value -1635.45
##
## F = -1635.45
## final value -1635.445748
## converged## Optimisation using starting value 8/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1459.9 |proj g|= 12
## At iterate 10 f = -1635.2 |proj g|= 4.8023
## At iterate 20 f = -1635.4 |proj g|= 0.053408
## At iterate 30 f = -1635.4 |proj g|= 0.26945
## At iterate 40 f = -1635.4 |proj g|= 0.083552
## At iterate 50 f = -1635.4 |proj g|= 0.015379
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1635.445724
## stopped after 52 iterations## Optimisation using starting value 9/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1370.5 |proj g|= 20
## At iterate 10 f = -1635.4 |proj g|= 1.5463
## At iterate 20 f = -1635.4 |proj g|= 0.069663
## At iterate 30 f = -1635.4 |proj g|= 0.29336
## At iterate 40 f = -1635.4 |proj g|= 0.16033
##
## iterations 48
## function evaluations 58
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0140004
## final function value -1635.45
##
## F = -1635.45
## final value -1635.445747
## converged## Optimisation using starting value 10/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1297.1 |proj g|= 20
## At iterate 10 f = -1635 |proj g|= 5.3865
## At iterate 20 f = -1635.4 |proj g|= 0.049719
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1635.441129
## stopped after 23 iterations## Optimisation using starting value 11/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1370.4 |proj g|= 20
## At iterate 10 f = -1635.4 |proj g|= 1.8538
##
## iterations 18
## function evaluations 37
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0164115
## final function value -1635.44
##
## F = -1635.44
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1635.441208
## converged## Optimisation using starting value 12/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1373.6 |proj g|= 20
## At iterate 10 f = -1635.4 |proj g|= 0.55454
##
## iterations 17
## function evaluations 23
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0255369
## final function value -1635.44
##
## F = -1635.44
## final value -1635.441240
## convergedprint(est.denver.model.B)## Optimisation for STmodel with 12 starting points.
## Results: 8 converged, 4 not converged, 0 failed.
## Best result for starting point 7, optimisation has converged
##
## No fixed parameters.
##
## Estimated parameters for all starting point(s):
## [,1] [,2] [,3]
## gamma.bc_st_no2 0.0062090995 0.006599608 0.006610820
## gamma.bc_st_smk 0.0522173639 0.051703033 0.051699867
## alpha.const.(Intercept) 0.1892293912 0.243936017 0.243776858
## alpha.const.tree_cover_100 0.0046762427 0.006591437 0.006562095
## alpha.const.impervious_2500 0.0174186970 0.013478089 0.013465414
## alpha.const.low_int_100 0.0120860605 0.014593970 0.014590474
## alpha.const.med_int_50 -0.0054002378 -0.008350045 -0.008333637
## alpha.const.high_int_100 0.0039325064 0.011433433 0.011432098
## alpha.const.pop_ct_1000 -0.0117617635 -0.010859412 -0.010844170
## alpha.const.dist_m_cafos -0.0218285758 -0.020202554 -0.020182808
## alpha.const.dist_m_og_wells -0.0027885270 0.006938867 0.006951907
## alpha.const.dist_m_wwtp 0.0017999181 -0.003781818 -0.003792484
## alpha.const.aadt_100 0.0347923990 0.020297694 0.020318340
## alpha.const.aadt_250 0.0223654732 0.022344864 0.022319732
## alpha.const.aadt_1000 -0.0072537280 -0.010120060 -0.010125223
## alpha.V1.(Intercept) 0.1642753582 0.178636679 0.178622552
## alpha.V1.tree_cover_100 -0.0129688194 -0.007844976 -0.007848563
## alpha.V1.impervious_2500 0.0130073263 0.003138609 0.003123109
## alpha.V1.low_int_100 -0.0040502127 -0.002775982 -0.002771018
## alpha.V1.med_int_50 -0.0134011261 -0.010220841 -0.010242218
## alpha.V1.high_int_100 0.0205863566 0.021306043 0.021314011
## alpha.V1.pop_ct_1000 -0.0102009204 -0.006645827 -0.006637243
## alpha.V1.dist_m_cafos 0.0031805666 -0.003804712 -0.003801116
## alpha.V1.dist_m_og_wells -0.0001515623 0.001987686 0.001972377
## alpha.V1.dist_m_wwtp -0.0017640557 -0.003374174 -0.003375790
## alpha.V1.aadt_100 -0.0252023998 -0.025189192 -0.025179605
## alpha.V1.aadt_250 -0.0230926273 -0.012513033 -0.012516054
## alpha.V1.aadt_1000 0.0029637176 -0.001501573 -0.001513485
## alpha.V2.(Intercept) 0.1614438120 0.183300320 0.183268769
## alpha.V2.tree_cover_100 -0.0039201823 -0.003182039 -0.003199660
## alpha.V2.impervious_2500 -0.0129678544 -0.014558645 -0.014580598
## alpha.V2.low_int_100 0.0034705004 0.004271867 0.004269628
## alpha.V2.med_int_50 -0.0088566551 -0.010601915 -0.010600557
## alpha.V2.high_int_100 -0.0108150580 -0.005522428 -0.005530511
## alpha.V2.pop_ct_1000 0.0082558630 0.008531245 0.008549593
## alpha.V2.dist_m_cafos 0.0100567351 0.010261381 0.010285709
## alpha.V2.dist_m_og_wells -0.0023454156 0.003669870 0.003656106
## alpha.V2.dist_m_wwtp -0.0004074542 -0.003195743 -0.003199130
## alpha.V2.aadt_100 -0.0196882673 -0.027581477 -0.027565941
## alpha.V2.aadt_250 0.0118562032 0.012766695 0.012762014
## alpha.V2.aadt_1000 -0.0047746198 -0.007479686 -0.007471430
## log.nugget.const.iid -15.0000000000 -12.301586463 -14.999996040
## log.nugget.V1.iid -14.6429024723 -7.106008778 -7.107788737
## log.nugget.V2.iid -15.0000000000 -7.620974222 -7.617071996
## nu.log.range.exp 0.0000000000 13.497616987 13.496689835
## nu.log.sill.exp -4.7059960711 -2.979560812 -2.979722590
## nu.log.nugget.(Intercept).exp -4.1963514833 -4.935190209 -4.934921194
## [,4] [,5] [,6]
## gamma.bc_st_no2 0.006606782 0.0057982734 0.006605847
## gamma.bc_st_smk 0.051703860 0.0529738986 0.051709804
## alpha.const.(Intercept) 0.243859851 0.2606675619 0.243908362
## alpha.const.tree_cover_100 0.006570133 0.0090419169 0.006563798
## alpha.const.impervious_2500 0.013471547 0.0150613829 0.013478401
## alpha.const.low_int_100 0.014591740 0.0156018309 0.014591659
## alpha.const.med_int_50 -0.008329657 -0.0101942492 -0.008340597
## alpha.const.high_int_100 0.011435480 0.0118084499 0.011432217
## alpha.const.pop_ct_1000 -0.010846383 -0.0122278851 -0.010838680
## alpha.const.dist_m_cafos -0.020181890 -0.0228376070 -0.020188625
## alpha.const.dist_m_og_wells 0.006946901 0.0055251356 0.006952237
## alpha.const.dist_m_wwtp -0.003785690 -0.0022115276 -0.003782381
## alpha.const.aadt_100 0.020317666 0.0174879994 0.020316336
## alpha.const.aadt_250 0.022323797 0.0243126042 0.022322717
## alpha.const.aadt_1000 -0.010124456 -0.0081782374 -0.010128323
## alpha.V1.(Intercept) 0.178647078 0.1785086700 0.178602916
## alpha.V1.tree_cover_100 -0.007846387 -0.0103300436 -0.007843408
## alpha.V1.impervious_2500 0.003128862 0.0039870707 0.003126096
## alpha.V1.low_int_100 -0.002770961 -0.0039995573 -0.002780179
## alpha.V1.med_int_50 -0.010242058 -0.0087999628 -0.010243653
## alpha.V1.high_int_100 0.021317005 0.0165687873 0.021309496
## alpha.V1.pop_ct_1000 -0.006637042 -0.0086001060 -0.006637788
## alpha.V1.dist_m_cafos -0.003796847 -0.0027752099 -0.003783661
## alpha.V1.dist_m_og_wells 0.001976021 0.0036886503 0.001961154
## alpha.V1.dist_m_wwtp -0.003375403 -0.0060100798 -0.003374056
## alpha.V1.aadt_100 -0.025184748 -0.0282904304 -0.025165225
## alpha.V1.aadt_250 -0.012519471 -0.0092431115 -0.012520104
## alpha.V1.aadt_1000 -0.001515458 0.0001764991 -0.001515384
## alpha.V2.(Intercept) 0.183287267 0.1918578438 0.183302278
## alpha.V2.tree_cover_100 -0.003194515 -0.0021046436 -0.003196501
## alpha.V2.impervious_2500 -0.014577253 -0.0121736760 -0.014573629
## alpha.V2.low_int_100 0.004269722 0.0045319340 0.004269193
## alpha.V2.med_int_50 -0.010597304 -0.0114719766 -0.010605240
## alpha.V2.high_int_100 -0.005528907 -0.0053398488 -0.005527698
## alpha.V2.pop_ct_1000 0.008549638 0.0064674459 0.008552768
## alpha.V2.dist_m_cafos 0.010285222 0.0084587404 0.010281139
## alpha.V2.dist_m_og_wells 0.003654501 0.0053575797 0.003655773
## alpha.V2.dist_m_wwtp -0.003196030 -0.0039596326 -0.003192532
## alpha.V2.aadt_100 -0.027568921 -0.0293141842 -0.027566782
## alpha.V2.aadt_250 0.012762727 0.0136781599 0.012760742
## alpha.V2.aadt_1000 -0.007470919 -0.0082915035 -0.007474092
## log.nugget.const.iid -15.000000000 -7.3607703383 -14.740913426
## log.nugget.V1.iid -7.109367887 -6.9920065975 -7.103233488
## log.nugget.V2.iid -7.618141184 -14.7044112710 -7.618123648
## nu.log.range.exp 13.499257735 13.3496770896 13.499896291
## nu.log.sill.exp -2.979349088 -2.9799999911 -2.979175398
## nu.log.nugget.(Intercept).exp -4.934814202 -4.9251779422 -4.934722135
## [,7] [,8] [,9]
## gamma.bc_st_no2 0.006608769 0.006611881 0.006606327
## gamma.bc_st_smk 0.051703365 0.051698402 0.051705019
## alpha.const.(Intercept) 0.243822987 0.243752207 0.243867384
## alpha.const.tree_cover_100 0.006566136 0.006560625 0.006572608
## alpha.const.impervious_2500 0.013471781 0.013463577 0.013475216
## alpha.const.low_int_100 0.014590699 0.014589960 0.014591235
## alpha.const.med_int_50 -0.008333643 -0.008333957 -0.008329567
## alpha.const.high_int_100 0.011432947 0.011431177 0.011435228
## alpha.const.pop_ct_1000 -0.010843227 -0.010844205 -0.010845792
## alpha.const.dist_m_cafos -0.020183805 -0.020182328 -0.020182016
## alpha.const.dist_m_og_wells 0.006949767 0.006952736 0.006945569
## alpha.const.dist_m_wwtp -0.003787490 -0.003794795 -0.003783885
## alpha.const.aadt_100 0.020317554 0.020318660 0.020317306
## alpha.const.aadt_250 0.022322542 0.022319003 0.022325757
## alpha.const.aadt_1000 -0.010126038 -0.010125320 -0.010125222
## alpha.V1.(Intercept) 0.178629397 0.178620464 0.178653318
## alpha.V1.tree_cover_100 -0.007844360 -0.007848343 -0.007841916
## alpha.V1.impervious_2500 0.003127615 0.003122332 0.003132623
## alpha.V1.low_int_100 -0.002772777 -0.002770030 -0.002770971
## alpha.V1.med_int_50 -0.010241952 -0.010241733 -0.010241023
## alpha.V1.high_int_100 0.021316517 0.021315011 0.021321299
## alpha.V1.pop_ct_1000 -0.006636606 -0.006636896 -0.006635892
## alpha.V1.dist_m_cafos -0.003794828 -0.003803356 -0.003794105
## alpha.V1.dist_m_og_wells 0.001971273 0.001972633 0.001976011
## alpha.V1.dist_m_wwtp -0.003372907 -0.003374837 -0.003371193
## alpha.V1.aadt_100 -0.025177600 -0.025179233 -0.025183228
## alpha.V1.aadt_250 -0.012520049 -0.012515929 -0.012523339
## alpha.V1.aadt_1000 -0.001515793 -0.001513267 -0.001517561
## alpha.V2.(Intercept) 0.183280585 0.183262792 0.183289758
## alpha.V2.tree_cover_100 -0.003196434 -0.003200743 -0.003192562
## alpha.V2.impervious_2500 -0.014576738 -0.014581488 -0.014574625
## alpha.V2.low_int_100 0.004269447 0.004269675 0.004269605
## alpha.V2.med_int_50 -0.010599889 -0.010600624 -0.010596310
## alpha.V2.high_int_100 -0.005529346 -0.005531116 -0.005528596
## alpha.V2.pop_ct_1000 0.008550674 0.008549425 0.008550333
## alpha.V2.dist_m_cafos 0.010284150 0.010285891 0.010284064
## alpha.V2.dist_m_og_wells 0.003654719 0.003656284 0.003653243
## alpha.V2.dist_m_wwtp -0.003195079 -0.003199793 -0.003192920
## alpha.V2.aadt_100 -0.027567470 -0.027565425 -0.027569892
## alpha.V2.aadt_250 0.012761913 0.012761953 0.012762642
## alpha.V2.aadt_1000 -0.007471315 -0.007471048 -0.007470133
## log.nugget.const.iid -15.000000000 -15.000000000 -15.000000000
## log.nugget.V1.iid -7.107260725 -7.107705894 -7.109151635
## log.nugget.V2.iid -7.617030091 -7.616346147 -7.617232539
## nu.log.range.exp 13.498729766 13.496128768 13.500518235
## nu.log.sill.exp -2.979686341 -2.979747813 -2.979529992
## nu.log.nugget.(Intercept).exp -4.934726995 -4.934852778 -4.934650953
## [,10] [,11] [,12]
## gamma.bc_st_no2 0.006596822 0.006593834 0.006593805
## gamma.bc_st_smk 0.051700890 0.051705743 0.051706151
## alpha.const.(Intercept) 0.243958173 0.244028589 0.244029652
## alpha.const.tree_cover_100 0.006601192 0.006605808 0.006606369
## alpha.const.impervious_2500 0.013477404 0.013484783 0.013485850
## alpha.const.low_int_100 0.014595348 0.014596198 0.014596162
## alpha.const.med_int_50 -0.008354777 -0.008354140 -0.008354013
## alpha.const.high_int_100 0.011434062 0.011435971 0.011436107
## alpha.const.pop_ct_1000 -0.010867322 -0.010866255 -0.010865858
## alpha.const.dist_m_cafos -0.020208824 -0.020209923 -0.020210132
## alpha.const.dist_m_og_wells 0.006934370 0.006931821 0.006931516
## alpha.const.dist_m_wwtp -0.003781403 -0.003774451 -0.003773551
## alpha.const.aadt_100 0.020290178 0.020289524 0.020289444
## alpha.const.aadt_250 0.022352934 0.022355654 0.022355994
## alpha.const.aadt_1000 -0.010116653 -0.010117277 -0.010117440
## alpha.V1.(Intercept) 0.178643850 0.178651584 0.178652489
## alpha.V1.tree_cover_100 -0.007847636 -0.007844738 -0.007844040
## alpha.V1.impervious_2500 0.003141753 0.003146001 0.003146798
## alpha.V1.low_int_100 -0.002775480 -0.002778207 -0.002778548
## alpha.V1.med_int_50 -0.010212518 -0.010213390 -0.010213515
## alpha.V1.high_int_100 0.021301469 0.021302089 0.021302507
## alpha.V1.pop_ct_1000 -0.006649812 -0.006649655 -0.006649516
## alpha.V1.dist_m_cafos -0.003812535 -0.003804350 -0.003803099
## alpha.V1.dist_m_og_wells 0.001996867 0.001995249 0.001994968
## alpha.V1.dist_m_wwtp -0.003376526 -0.003375663 -0.003375195
## alpha.V1.aadt_100 -0.025197894 -0.025196505 -0.025196218
## alpha.V1.aadt_250 -0.012508341 -0.012511674 -0.012512344
## alpha.V1.aadt_1000 -0.001495058 -0.001497345 -0.001497910
## alpha.V2.(Intercept) 0.183302241 0.183319879 0.183320643
## alpha.V2.tree_cover_100 -0.003177040 -0.003173366 -0.003172861
## alpha.V2.impervious_2500 -0.014553306 -0.014549452 -0.014548807
## alpha.V2.low_int_100 0.004272952 0.004272712 0.004272619
## alpha.V2.med_int_50 -0.010601788 -0.010601226 -0.010601091
## alpha.V2.high_int_100 -0.005520156 -0.005518543 -0.005518373
## alpha.V2.pop_ct_1000 0.008522607 0.008524045 0.008524342
## alpha.V2.dist_m_cafos 0.010253588 0.010252525 0.010252432
## alpha.V2.dist_m_og_wells 0.003676225 0.003674680 0.003674291
## alpha.V2.dist_m_wwtp -0.003198048 -0.003194030 -0.003193257
## alpha.V2.aadt_100 -0.027586892 -0.027588508 -0.027588669
## alpha.V2.aadt_250 0.012768977 0.012768874 0.012768822
## alpha.V2.aadt_1000 -0.007482735 -0.007483052 -0.007483005
## log.nugget.const.iid -11.986919760 -11.999339558 -12.003766248
## log.nugget.V1.iid -7.106592398 -7.106151711 -7.106031383
## log.nugget.V2.iid -7.622729793 -7.623526874 -7.623464700
## nu.log.range.exp 13.496344562 13.498723296 13.499031877
## nu.log.sill.exp -2.979484160 -2.979266929 -2.979294251
## nu.log.nugget.(Intercept).exp -4.935379234 -4.935196299 -4.935104756
##
## Function value(s):
## [1] 1393.212 1635.442 1635.446 1635.446 1629.644 1635.446 1635.446 1635.446
## [9] 1635.446 1635.441 1635.441 1635.4413.2.3 Cross-validation
Define the CV groups
set.seed(1000)
site_idsB <- colnames(denver.model.B$D.beta)[which(str_detect(colnames(denver.model.B$D.beta), "d_") == T)]
site_idsB## [1] "d_1" "d_10" "d_11" "d_12" "d_13" "d_14" "d_15" "d_16" "d_17" "d_18"
## [11] "d_19" "d_2" "d_20" "d_21" "d_22" "d_23" "d_25" "d_26" "d_28" "d_29"
## [21] "d_3" "d_31" "d_32" "d_33" "d_34" "d_35" "d_36" "d_37" "d_38" "d_39"
## [31] "d_4" "d_40" "d_41" "d_42" "d_43" "d_44" "d_45" "d_46" "d_47" "d_48"
## [41] "d_49" "d_5" "d_50" "d_51" "d_52" "d_53" "d_54" "d_55" "d_56" "d_57"
## [51] "d_58" "d_6" "d_60" "d_61" "d_62" "d_63" "d_64" "d_65" "d_66" "d_67"
## [61] "d_68" "d_7" "d_8"Ind.cv.B <- createCV(denver.model.B, groups = 10, #min.dist = .1,
subset = site_idsB)
ID.cv.B <- sapply(split(denver.model.B$obs$ID, Ind.cv.B), unique)
print(sapply(ID.cv.B, length))## 0 1 2 3 4 5 6 7 8 9 10
## 1 7 7 7 6 6 6 6 6 6 6table(Ind.cv.B)## Ind.cv.B
## 0 1 2 3 4 5 6 7 8 9 10
## 231 75 80 97 73 64 93 80 81 66 81I.col.B <- apply(sapply(ID.cv.B, function(x) denver.model.B$locations$ID%in% x), 1,
function(x) if(sum(x)==1) which(x) else 0)
names(I.col.B) <- denver.model.B$locations$ID
print(I.col.B)## central d_1 d_10 d_11 d_12 d_13 d_14 d_15 d_16 d_17
## 1 6 2 3 4 4 4 9 4 10
## d_18 d_19 d_2 d_20 d_21 d_22 d_23 d_25 d_26 d_28
## 8 5 11 3 2 11 4 9 11 9
## d_29 d_3 d_31 d_32 d_33 d_34 d_35 d_36 d_37 d_38
## 2 4 6 8 5 5 10 8 3 8
## d_39 d_4 d_40 d_41 d_42 d_43 d_44 d_45 d_46 d_47
## 6 7 7 3 7 8 11 3 9 5
## d_48 d_49 d_5 d_50 d_51 d_52 d_53 d_54 d_55 d_56
## 9 4 6 9 10 3 7 5 7 6
## d_57 d_58 d_6 d_60 d_61 d_62 d_63 d_64 d_65 d_66
## 2 2 11 5 2 11 10 3 8 6
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 7 10 2 10 0 0 0 0 0par(mfrow=c(1,1))
plot(denver.model.B$locations$long,
denver.model.B$locations$lat,
pch=23+floor(I.col.B/max(I.col.B)+.5), bg=I.col.B,
xlab="Longitude", ylab="Latitude")
map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)## [[1]]
## NULLID starting conditions using previous model without CV:
x.init.B.cv <- coef(est.denver.model.B, pars="cov")[,c("par","init")]
x.init.B.cvRun the model with cross validation
est.denver.B.cv <- estimateCV(denver.model.B, x.init.B.cv, Ind.cv.B)##
## ***************************## Estimation of CV-set 1/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1485.5 |proj g|= 12.41
## At iterate 10 f = -1486.4 |proj g|= 0.11533
##
## iterations 15
## function evaluations 25
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0355727
## final function value -1486.45
##
## F = -1486.45
## final value -1486.450249
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1158.2 |proj g|= 20
## At iterate 10 f = -1483.6 |proj g|= 7.9761
## At iterate 20 f = -1486.4 |proj g|= 1.0742
## At iterate 30 f = -1486.4 |proj g|= 1.1381
## At iterate 40 f = -1486.4 |proj g|= 0.41022
## At iterate 50 f = -1486.4 |proj g|= 0.03182
##
## iterations 52
## function evaluations 62
## segments explored during Cauchy searches 57
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0186954
## final function value -1486.45
##
## F = -1486.45
## final value -1486.449712
## converged## ##
## ***************************## Estimation of CV-set 2/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1499.5 |proj g|= 8.6836
## At iterate 10 f = -1499.9 |proj g|= 0.34583
##
## iterations 16
## function evaluations 40
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0256119
## final function value -1499.86
##
## F = -1499.86
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1499.856343
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1176.8 |proj g|= 20
## At iterate 10 f = -1498.3 |proj g|= 4.3979
## At iterate 20 f = -1499.7 |proj g|= 0.16062
## At iterate 30 f = -1499.8 |proj g|= 1.9279
## At iterate 40 f = -1499.8 |proj g|= 0.51062
## At iterate 50 f = -1499.9 |proj g|= 0.050572
## At iterate 60 f = -1499.9 |proj g|= 0.03524
##
## iterations 66
## function evaluations 81
## segments explored during Cauchy searches 71
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0268738
## final function value -1499.86
##
## F = -1499.86
## final value -1499.856374
## converged## ##
## ***************************## Estimation of CV-set 3/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1468.7 |proj g|= 10.065
## At iterate 10 f = -1469.1 |proj g|= 0.1826
##
## iterations 12
## function evaluations 16
## segments explored during Cauchy searches 13
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0485505
## final function value -1469.05
##
## F = -1469.05
## final value -1469.051766
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1150.3 |proj g|= 20
## At iterate 10 f = -1466.2 |proj g|= 8.1879
## At iterate 20 f = -1469 |proj g|= 0.52463
## At iterate 30 f = -1469 |proj g|= 0.83975
## At iterate 40 f = -1469.1 |proj g|= 0.072265
##
## iterations 45
## function evaluations 55
## segments explored during Cauchy searches 50
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00180551
## final function value -1469.05
##
## F = -1469.05
## final value -1469.051878
## converged## ##
## ***************************## Estimation of CV-set 4/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1481.8 |proj g|= 17.401
## At iterate 10 f = -1482.6 |proj g|= 0.10468
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1482.621115
## stopped after 12 iterations## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1162.2 |proj g|= 20
## At iterate 10 f = -1478.9 |proj g|= 8.2947
## At iterate 20 f = -1482.6 |proj g|= 0.17772
## At iterate 30 f = -1482.6 |proj g|= 0.15057
## At iterate 40 f = -1482.6 |proj g|= 0.17725
## At iterate 50 f = -1482.6 |proj g|= 0.037305
##
## iterations 51
## function evaluations 64
## segments explored during Cauchy searches 56
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0317541
## final function value -1482.62
##
## F = -1482.62
## final value -1482.618825
## converged## ##
## ***************************## Estimation of CV-set 5/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1505.6 |proj g|= 11.049
## At iterate 10 f = -1506.1 |proj g|= 0.0311
##
## iterations 11
## function evaluations 16
## segments explored during Cauchy searches 11
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0247902
## final function value -1506.14
##
## F = -1506.14
## final value -1506.140050
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1186 |proj g|= 20
## At iterate 10 f = -1505.4 |proj g|= 3.5786
## At iterate 20 f = -1506 |proj g|= 0.20557
## At iterate 30 f = -1506.1 |proj g|= 0.51914
## At iterate 40 f = -1506.1 |proj g|= 0.04747
##
## iterations 42
## function evaluations 62
## segments explored during Cauchy searches 47
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0181309
## final function value -1506.14
##
## F = -1506.14
## final value -1506.140042
## converged## ##
## ***************************## Estimation of CV-set 6/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1478.5 |proj g|= 10.065
## At iterate 10 f = -1478.9 |proj g|= 0.34451
##
## iterations 14
## function evaluations 18
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0260103
## final function value -1478.86
##
## F = -1478.86
## final value -1478.857207
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1146.8 |proj g|= 20
## At iterate 10 f = -1477.8 |proj g|= 3.6991
## At iterate 20 f = -1478.8 |proj g|= 0.15713
## At iterate 30 f = -1478.8 |proj g|= 0.17528
## At iterate 40 f = -1478.9 |proj g|= 0.080178
## At iterate 50 f = -1478.9 |proj g|= 0.11838
##
## iterations 55
## function evaluations 74
## segments explored during Cauchy searches 60
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0120084
## final function value -1478.86
##
## F = -1478.86
## final value -1478.857194
## converged## ##
## ***************************## Estimation of CV-set 7/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1491.1 |proj g|= 2.2904
## At iterate 10 f = -1492.3 |proj g|= 0.19837
##
## iterations 16
## function evaluations 20
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0734352
## final function value -1492.31
##
## F = -1492.31
## final value -1492.306722
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1156.4 |proj g|= 20
## At iterate 10 f = -1490.3 |proj g|= 2.4523
## At iterate 20 f = -1492 |proj g|= 0.42866
## At iterate 30 f = -1492.3 |proj g|= 0.62506
## At iterate 40 f = -1492.3 |proj g|= 0.058579
##
## iterations 46
## function evaluations 66
## segments explored during Cauchy searches 51
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.101364
## final function value -1492.31
##
## F = -1492.31
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1492.306663
## converged## ##
## ***************************## Estimation of CV-set 8/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1478.8 |proj g|= 9.4926
## At iterate 10 f = -1479.8 |proj g|= 1.0155
##
## iterations 15
## function evaluations 17
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0169185
## final function value -1479.84
##
## F = -1479.84
## final value -1479.840150
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1151.7 |proj g|= 20
## At iterate 10 f = -1478.7 |proj g|= 2.4318
## At iterate 20 f = -1479.8 |proj g|= 0.48086
## At iterate 30 f = -1479.9 |proj g|= 0.028233
## At iterate 40 f = -1479.9 |proj g|= 0.23195
## At iterate 50 f = -1479.9 |proj g|= 0.0444
##
## iterations 50
## function evaluations 68
## segments explored during Cauchy searches 55
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0444004
## final function value -1479.88
##
## F = -1479.88
## final value -1479.880325
## converged## ##
## ***************************## Estimation of CV-set 9/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1501.1 |proj g|= 11.422
## At iterate 10 f = -1501.6 |proj g|= 0.10687
##
## iterations 17
## function evaluations 38
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0211563
## final function value -1501.55
##
## F = -1501.55
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1501.552882
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1174.8 |proj g|= 20
## At iterate 10 f = -1499.8 |proj g|= 4.3954
## At iterate 20 f = -1501.5 |proj g|= 0.10567
## At iterate 30 f = -1501.5 |proj g|= 0.31642
## At iterate 40 f = -1501.5 |proj g|= 0.056127
## At iterate 50 f = -1501.6 |proj g|= 0.026776
##
## iterations 51
## function evaluations 68
## segments explored during Cauchy searches 56
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0267776
## final function value -1501.55
##
## F = -1501.55
## final value -1501.552897
## converged## ##
## ***************************## Estimation of CV-set 10/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 1 variables are exactly at the bounds
## At iterate 0 f= -1497.5 |proj g|= 10.065
## At iterate 10 f = -1498.7 |proj g|= 0.11568
##
## iterations 15
## function evaluations 27
## segments explored during Cauchy searches 16
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.017471
## final function value -1498.69
##
## F = -1498.69
## final value -1498.694111
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1162.3 |proj g|= 20
## At iterate 10 f = -1497.8 |proj g|= 2.9104
## At iterate 20 f = -1498.7 |proj g|= 0.23363
## At iterate 30 f = -1498.7 |proj g|= 0.26079
## At iterate 40 f = -1498.7 |proj g|= 0.089778
## At iterate 50 f = -1498.7 |proj g|= 0.044009
## At iterate 60 f = -1498.7 |proj g|= 0.066407
##
## iterations 60
## function evaluations 83
## segments explored during Cauchy searches 65
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0664066
## final function value -1498.69
##
## F = -1498.69
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1498.694142
## converged## print(est.denver.B.cv)## Cross-validation parameter estimation for STmodel
## with 10 CV-groups and 2 starting points.
## Results: 6 converged, 4 not converged.
##
## No fixed parameters.
##
## Estimated function values and convergence info:
## value convergence conv eigen.min eigen.all.min
## 1 1486.450 TRUE FALSE -0.00105621694 NA
## 2 1499.856 TRUE TRUE 0.00249978947 NA
## 3 1469.052 TRUE FALSE -0.00185320296 NA
## 4 1482.619 TRUE TRUE 0.00247626817 NA
## 5 1506.140 TRUE TRUE 0.00046616451 NA
## 6 1478.857 TRUE TRUE 0.00003777341 NA
## 7 1492.307 TRUE TRUE 0.00344725005 NA
## 8 1479.880 TRUE TRUE 0.08593726122 NA
## 9 1501.553 TRUE FALSE -0.00104913011 NA
## 10 1498.694 TRUE FALSE -0.00029446967 NApar(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.B, pars="all"),
plotCI((1:length(par))+.3, par, uiw=1.96*sd,
col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.B.cv, "all", boxwex=.4, col="grey", add=TRUE)
#' Save the results as an .rdata object
save(denver.data.B, denver.model.B, est.denver.model.B, est.denver.B.cv,
file = here::here("Results", "Denver_ST_Model_B.rdata"))3.2.4 Prediction using the CV model
Making predictions using the CV model. Printing out the CV summary statistics as well
pred.B.cv <- predictCV(denver.model.B, est.denver.B.cv, LTA = T)
pred.B.cv.log <- predictCV(denver.model.B, est.denver.B.cv,
LTA = T, transform="unbiased")
head(pred.B.cv$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 0.2210964 0.3407053 0.4451902 0.3810322 0.3209532 0.3787595
## 2009-01-05 NA 0.2361515 0.3502648 0.4450109 0.3852111 0.3217106 0.3741248
## 2009-01-12 NA 0.2940590 0.3976192 0.4786138 0.4148080 0.3480191 0.3948416
## 2009-01-19 NA 0.2467917 0.3516863 0.4297334 0.3799835 0.3101962 0.3510143
## 2009-01-26 NA 0.2461414 0.3466328 0.4178441 0.3718647 0.2995256 0.3337137
## 2009-02-02 NA 0.2552510 0.3497491 0.4140709 0.3677949 0.2934817 0.3202370
## d_15 d_16 d_17 d_18 d_19 d_2
## 2008-12-29 0.1843625 0.4566230 0.2142486 0.3912299 0.2605300 0.3811699
## 2009-01-05 0.2039089 0.4686705 0.2351228 0.3928683 0.2811398 0.4117604
## 2009-01-12 0.2596049 0.5067382 0.2920126 0.4311673 0.3468917 0.4869437
## 2009-01-19 0.2242828 0.4816472 0.2582672 0.3792863 0.3020073 0.4559239
## 2009-01-26 0.2271822 0.4851848 0.2626188 0.3672268 0.3062898 0.4748955
## 2009-02-02 0.2363035 0.4952447 0.2732028 0.3633895 0.3209703 0.5067628
## d_20 d_21 d_22 d_23 d_25 d_26
## 2008-12-29 0.4234567 0.2606050 0.2275603 0.4349549 0.4055532 0.3411497
## 2009-01-05 0.4240978 0.2733945 0.2471630 0.4409138 0.4040741 0.3500121
## 2009-01-12 0.4583904 0.3239300 0.3105965 0.4725566 0.4392007 0.4028886
## 2009-01-19 0.4099209 0.2810676 0.2662320 0.4403390 0.3843216 0.3483784
## 2009-01-26 0.3980150 0.2789123 0.2694307 0.4356831 0.3692330 0.3420712
## 2009-02-02 0.3936739 0.2847026 0.2824028 0.4361784 0.3624013 0.3463683
## d_28 d_29 d_3 d_31 d_32 d_33
## 2008-12-29 0.2566995 0.1907858 0.3518235 0.5065135 0.1853651 0.1994813
## 2009-01-05 0.2645547 0.2024712 0.3592448 0.5056800 0.1969113 0.2138933
## 2009-01-12 0.3091322 0.2515903 0.3920378 0.5482645 0.2449154 0.2735540
## 2009-01-19 0.2639321 0.2066502 0.3603044 0.4869079 0.2022887 0.2228176
## 2009-01-26 0.2588621 0.2014073 0.3551129 0.4740706 0.1987845 0.2216202
## 2009-02-02 0.2624807 0.2028091 0.3537572 0.4734546 0.2025953 0.2313039
## d_34 d_35 d_36 d_37 d_38 d_39
## 2008-12-29 0.3969433 0.3569734 0.4090694 0.3445539 0.3249153 0.4141556
## 2009-01-05 0.4040162 0.3604100 0.4078161 0.3472050 0.3351181 0.4130795
## 2009-01-12 0.4563984 0.3999580 0.4433621 0.3835707 0.3819081 0.4550802
## 2009-01-19 0.3985276 0.3491112 0.3890289 0.3373043 0.3383443 0.3924162
## 2009-01-26 0.3904233 0.3367458 0.3749800 0.3277981 0.3343279 0.3771649
## 2009-02-02 0.3934978 0.3311190 0.3697514 0.3261090 0.3381751 0.3727092
## d_4 d_40 d_41 d_42 d_43 d_44
## 2008-12-29 0.2325459 0.2325977 0.4548838 0.2489322 0.3746541 0.2797726
## 2009-01-05 0.2590560 0.2532798 0.4463912 0.2706392 0.3749866 0.2914816
## 2009-01-12 0.3613801 0.3495479 0.4716445 0.3681047 0.4119964 0.3471394
## 2009-01-19 0.2828344 0.2644740 0.4143547 0.2845927 0.3588645 0.2952662
## 2009-01-26 0.2892164 0.2636114 0.3939727 0.2858482 0.3456141 0.2913724
## 2009-02-02 0.3173627 0.2835933 0.3816048 0.3086646 0.3406637 0.2977930
## d_45 d_46 d_47 d_48 d_49 d_5
## 2008-12-29 0.4413909 0.2149524 0.2609522 0.4122577 0.3977901 0.3342450
## 2009-01-05 0.4348121 0.2284492 0.2769073 0.4133927 0.3959995 0.3444773
## 2009-01-12 0.4617900 0.2786291 0.3380616 0.4511285 0.4197358 0.3975420
## 2009-01-19 0.4058195 0.2389364 0.2887105 0.3988431 0.3792947 0.3453992
## 2009-01-26 0.3861350 0.2392230 0.2887312 0.3863215 0.3659383 0.3398281
## 2009-02-02 0.3736618 0.2480004 0.2994167 0.3820209 0.3571229 0.3439617
## d_50 d_51 d_52 d_53 d_54 d_55
## 2008-12-29 0.2218148 0.3477089 0.4647677 0.1733572 0.2977684 0.2618503
## 2009-01-05 0.2325491 0.3496273 0.4611208 0.2072953 0.3092335 0.2782711
## 2009-01-12 0.2800061 0.3877569 0.4911496 0.3170293 0.3659009 0.3703172
## 2009-01-19 0.2376805 0.3357063 0.4384762 0.2458396 0.3120790 0.2811137
## 2009-01-26 0.2354743 0.3224670 0.4224626 0.2594877 0.3076579 0.2762669
## 2009-02-02 0.2419415 0.3163923 0.4141403 0.2947798 0.3139416 0.2924548
## d_56 d_57 d_58 d_6 d_60 d_61
## 2008-12-29 0.2990110 0.3838416 0.2838885 0.3271021 0.3373696 0.3233478
## 2009-01-05 0.3111394 0.3819728 0.2974450 0.3331940 0.3456672 0.3352979
## 2009-01-12 0.3661266 0.4178219 0.3486795 0.3833369 0.3991938 0.3849843
## 2009-01-19 0.3159585 0.3602437 0.3063699 0.3261771 0.3422947 0.3412537
## 2009-01-26 0.3124399 0.3433422 0.3045432 0.3173514 0.3348964 0.3382019
## 2009-02-02 0.3187252 0.3343562 0.3103726 0.3193008 0.3383342 0.3430597
## d_62 d_63 d_64 d_65 d_66 d_67
## 2008-12-29 0.2519884 0.3333118 0.4327805 0.2995690 0.2630535 0.2956804
## 2009-01-05 0.2683571 0.3393123 0.4244691 0.3106226 0.2755992 0.3163226
## 2009-01-12 0.3286029 0.3816060 0.4496926 0.3581327 0.3311175 0.4130007
## 2009-01-19 0.2811554 0.3338855 0.3919249 0.3150101 0.2817213 0.3292919
## 2009-01-26 0.2814344 0.3252260 0.3703794 0.3110076 0.2793423 0.3312551
## 2009-02-02 0.2916994 0.3240504 0.3559642 0.3143172 0.2872409 0.3559454
## d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.2503416 0.3365727 0.1628510 NA NA NA NA NA
## 2009-01-05 0.2656210 0.3479577 0.1836460 NA NA NA NA NA
## 2009-01-12 0.3171798 0.3970807 0.2407996 NA NA NA NA NA
## 2009-01-19 0.2786753 0.3527917 0.2080467 NA NA NA NA NA
## 2009-01-26 0.2791480 0.3491896 0.2145052 NA NA NA NA NA
## 2009-02-02 0.2869906 0.3535080 0.2286322 NA NA NA NA NAtail(pred.B.cv$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 1.0113268 1.0829636 1.0081086 1.0859559 0.9650805 0.9735064
## 2020-11-09 NA 0.6223434 0.6980942 0.6201038 0.6950385 0.5700992 0.6048310
## 2020-11-16 NA 1.0213745 1.1103616 1.0426701 1.1075152 0.9910148 1.0294706
## 2020-11-23 NA 0.9513426 1.0393384 0.9811967 1.0335056 0.9210381 0.9698794
## 2020-11-30 NA 1.1462592 1.2409710 1.1980696 1.2411175 1.1358406 1.1928841
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2
## 2020-11-02 1.1155383 1.1519096 1.1656897 1.0205160 1.1198417 1.3074259
## 2020-11-09 0.7215591 0.7274304 0.7610358 0.6426825 0.7283319 0.8470043
## 2020-11-16 1.1226782 1.1281099 1.1560897 1.0668613 1.1272710 1.2226713
## 2020-11-23 1.0365859 1.0187177 1.0711838 1.0075477 1.0444774 1.1106700
## 2020-11-30 1.2228471 1.2138336 1.2602210 1.2234996 1.2272958 1.2697088
## 2020-12-07 NA NA NA NA NA NA
## d_20 d_21 d_22 d_23 d_25 d_26
## 2020-11-02 1.0164026 1.0525841 1.0963512 1.1078621 0.9961886 1.0568900
## 2020-11-09 0.6374244 0.6597018 0.6943842 0.7092348 0.6272813 0.6582018
## 2020-11-16 1.0655825 1.0656339 1.0938027 1.1214530 1.0575167 1.0817807
## 2020-11-23 1.0023643 0.9920130 1.0160235 1.0352619 1.0012920 1.0204201
## 2020-11-30 1.2247127 1.1901809 1.2076444 1.2335757 1.2228911 1.2327153
## 2020-12-07 NA NA NA NA NA NA
## d_28 d_29 d_3 d_31 d_32 d_33
## 2020-11-02 0.9279066 1.0292349 1.1116712 0.9986040 0.9765376 0.9604948
## 2020-11-09 0.5427529 0.6485138 0.7118513 0.6279082 0.5851825 0.5736401
## 2020-11-16 0.9508429 1.0584710 1.1226925 1.0558598 0.9949439 0.9757666
## 2020-11-23 0.8770509 0.9913678 1.0429693 0.9979033 0.9206461 0.9060979
## 2020-11-30 1.0790052 1.1963793 1.2468834 1.2098419 1.1227990 1.1019361
## 2020-12-07 NA NA NA NA NA NA
## d_34 d_35 d_36 d_37 d_38 d_39
## 2020-11-02 1.0681293 1.0803513 0.9744032 0.9643484 1.0729829 0.9583423
## 2020-11-09 0.6897761 0.7054644 0.6003245 0.5786706 0.6812147 0.6021643
## 2020-11-16 1.1099127 1.1320171 1.0275643 0.9987399 1.0922749 1.0233930
## 2020-11-23 1.0497390 1.0743833 0.9704524 0.9298452 1.0167910 0.9766798
## 2020-11-30 1.2610199 1.2945195 1.1887477 1.1399516 1.2206827 1.1933249
## 2020-12-07 NA NA NA NA NA NA
## d_4 d_40 d_41 d_42 d_43 d_44
## 2020-11-02 1.0371786 0.8953660 0.9363574 1.0151816 0.9861491 1.0307745
## 2020-11-09 0.6345221 0.5372357 0.5697597 0.6275572 0.6119862 0.6397091
## 2020-11-16 1.0422628 0.9309535 1.0106574 1.0458053 1.0392250 1.0525235
## 2020-11-23 0.9827068 0.8692533 0.9586908 0.9926211 0.9793235 0.9866129
## 2020-11-30 1.1758109 1.0646080 1.1904373 1.1924960 1.1983007 1.1905183
## 2020-12-07 NA NA NA NA NA NA
## d_45 d_46 d_47 d_48 d_49 d_5
## 2020-11-02 0.9816400 0.9600055 1.0517271 1.0304915 1.0142244 1.0429725
## 2020-11-09 0.6166026 0.5626067 0.6692316 0.6569391 0.6306741 0.6677092
## 2020-11-16 1.0546074 0.9594819 1.0736460 1.0792439 1.0574433 1.0689086
## 2020-11-23 1.0044121 0.8804099 1.0005365 1.0221229 0.9897138 1.0061854
## 2020-11-30 1.2341302 1.0725935 1.1852746 1.2374049 1.2047201 1.2084411
## 2020-12-07 NA NA NA NA NA NA
## d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 0.9220988 1.0251997 1.0167359 1.1627859 1.0316722 0.9869116
## 2020-11-09 0.5330091 0.6585428 0.6399447 0.7298017 0.6503707 0.6183035
## 2020-11-16 0.9321225 1.0861155 1.0712913 1.1307733 1.0658968 1.0454210
## 2020-11-23 0.8568012 1.0295448 1.0166389 1.0568680 1.0014371 1.0014890
## 2020-11-30 1.0515841 1.2420104 1.2395524 1.2345028 1.2108328 1.2113879
## 2020-12-07 NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 1.0384630 0.9957854 1.1107714 0.9927483 1.024133 1.088869 1.0521552
## 2020-11-09 0.6551031 0.6371126 0.7153502 0.6114409 0.646523 0.698486 0.6584805
## 2020-11-16 1.0558113 1.0647092 1.1226133 1.0350975 1.061061 1.100196 1.0563271
## 2020-11-23 0.9915705 1.0173536 1.0486268 0.9746616 1.000040 1.030916 0.9854171
## 2020-11-30 1.1917756 1.2401399 1.2456075 1.1860702 1.208811 1.231508 1.1811265
## 2020-12-07 NA NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68
## 2020-11-02 1.0454855 0.9504433 1.0868516 1.0057229 0.9949267 1.0997216
## 2020-11-09 0.6644575 0.5920252 0.6806218 0.6155059 0.6127107 0.6954856
## 2020-11-16 1.0808183 1.0294935 1.0958678 1.0237729 1.0328983 1.0984914
## 2020-11-23 1.0175663 0.9847705 1.0278611 0.9507318 0.9723737 1.0191424
## 2020-11-30 1.2273449 1.2160036 1.2377496 1.1513341 1.1616442 1.2126704
## 2020-12-07 NA NA NA NA NA NA
## d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.0963415 1.0403187 NA NA NA NA NA
## 2020-11-09 0.7056589 0.6301959 NA NA NA NA NA
## 2020-11-16 1.1100093 1.0197611 NA NA NA NA NA
## 2020-11-23 1.0400124 0.9290728 NA NA NA NA NA
## 2020-11-30 1.2421562 1.1119611 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA NAhead(pred.B.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 1.280886 1.442882 1.603164 1.504168 1.416460 1.500753
## 2009-01-05 NA 1.300018 1.456520 1.602472 1.510060 1.417152 1.493412
## 2009-01-12 NA 1.377243 1.526948 1.656867 1.555054 1.454587 1.524313
## 2009-01-19 NA 1.313434 1.458230 1.577528 1.501533 1.400318 1.458659
## 2009-01-26 NA 1.312396 1.450742 1.558642 1.489152 1.385232 1.433410
## 2009-02-02 NA 1.324263 1.455163 1.552587 1.482922 1.376717 1.414048
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2008-12-29 1.234092 1.622277 1.272186 1.518502 1.331824 1.504657 1.568698
## 2009-01-05 1.258280 1.641498 1.298708 1.520634 1.359256 1.550922 1.569307
## 2009-01-12 1.330190 1.704788 1.374442 1.579676 1.451349 1.671572 1.623696
## 2009-01-19 1.283894 1.662214 1.328597 1.499549 1.387416 1.620149 1.546581
## 2009-01-26 1.287512 1.667836 1.334197 1.481366 1.393180 1.650878 1.528040
## 2009-02-02 1.299223 1.684493 1.348244 1.475535 1.413634 1.704099 1.521239
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2008-12-29 1.331814 1.290404 1.587503 1.539605 1.445630 1.326670 1.242000
## 2009-01-05 1.348752 1.315547 1.596562 1.537121 1.458052 1.336951 1.256408
## 2009-01-12 1.418474 1.401324 1.647500 1.591883 1.536811 1.397730 1.319485
## 2009-01-19 1.358804 1.340210 1.594950 1.506721 1.454952 1.335822 1.261356
## 2009-01-26 1.355750 1.344259 1.587286 1.484032 1.445540 1.328953 1.254641
## 2009-02-02 1.363522 1.361621 1.587877 1.473829 1.451564 1.333682 1.256308
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2008-12-29 1.460868 1.703985 1.235975 1.252950 1.526477 1.467355 1.545835
## 2009-01-05 1.471354 1.702174 1.250033 1.270856 1.536971 1.472051 1.543535
## 2009-01-12 1.520046 1.775867 1.311233 1.348720 1.619311 1.531110 1.599058
## 2009-01-19 1.472273 1.669895 1.256295 1.281784 1.528005 1.454944 1.514230
## 2009-01-26 1.464414 1.648365 1.251725 1.280076 1.515465 1.436855 1.492896
## 2009-02-02 1.462251 1.647172 1.256370 1.292396 1.519971 1.428635 1.484952
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2008-12-29 1.449681 1.421070 1.553657 1.295786 1.295853 1.618781 1.317194
## 2009-01-05 1.453162 1.435305 1.551630 1.330206 1.322545 1.604686 1.345704
## 2009-01-12 1.506645 1.503748 1.617861 1.473146 1.455818 1.645360 1.483086
## 2009-01-19 1.438254 1.439394 1.519330 1.361568 1.336797 1.553453 1.363964
## 2009-01-26 1.424426 1.433424 1.496124 1.370043 1.335408 1.521876 1.365436
## 2009-02-02 1.421852 1.438796 1.489312 1.408963 1.362177 1.502990 1.396760
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2008-12-29 1.493539 1.359569 1.597085 1.272426 1.332386 1.549962 1.529587
## 2009-01-05 1.493684 1.375161 1.586212 1.289541 1.353515 1.551512 1.526440
## 2009-01-12 1.549681 1.453479 1.629225 1.355738 1.438590 1.610985 1.562736
## 2009-01-19 1.469236 1.379692 1.540250 1.302846 1.369090 1.528761 1.500499
## 2009-01-26 1.449693 1.374080 1.509994 1.303108 1.368932 1.509610 1.480353
## 2009-02-02 1.442381 1.382739 1.491099 1.314509 1.383491 1.503031 1.467181
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2008-12-29 1.434334 1.281188 1.453823 1.634860 1.221315 1.382354 1.334320
## 2009-01-05 1.448754 1.294839 1.456264 1.628497 1.263105 1.397983 1.356013
## 2009-01-12 1.527399 1.357607 1.512543 1.677768 1.409238 1.479202 1.486371
## 2009-01-19 1.449549 1.301210 1.435570 1.591380 1.312117 1.401460 1.359227
## 2009-01-26 1.441294 1.298233 1.416484 1.565857 1.329913 1.395088 1.352416
## 2009-02-02 1.447108 1.306569 1.407750 1.552694 1.377501 1.403733 1.374302
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2008-12-29 1.384677 1.506485 1.363188 1.425464 1.438195 1.418053 1.322314
## 2009-01-05 1.401252 1.503444 1.381584 1.433735 1.449856 1.434883 1.343726
## 2009-01-12 1.480161 1.558110 1.454018 1.507055 1.529278 1.507776 1.426785
## 2009-01-19 1.407495 1.470762 1.393623 1.423006 1.444452 1.443096 1.360360
## 2009-01-26 1.402355 1.445976 1.390948 1.410245 1.433610 1.438563 1.360492
## 2009-02-02 1.411045 1.432935 1.398977 1.412801 1.438395 1.445461 1.374339
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2008-12-29 1.433043 1.583393 1.385504 1.335772 1.380232 1.318942 1.436932
## 2009-01-05 1.441320 1.569890 1.400573 1.352325 1.408606 1.338926 1.453164
## 2009-01-12 1.503268 1.609634 1.468417 1.429239 1.551188 1.409472 1.526126
## 2009-01-19 1.432959 1.518997 1.406196 1.360122 1.426315 1.355989 1.459843
## 2009-01-26 1.420398 1.486390 1.400383 1.356700 1.428866 1.356434 1.454456
## 2009-02-02 1.418572 1.464942 1.404875 1.367311 1.464387 1.366962 1.460643
## d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 1.208451 NA NA NA NA NA
## 2009-01-05 1.233546 NA NA NA NA NA
## 2009-01-12 1.305825 NA NA NA NA NA
## 2009-01-19 1.263522 NA NA NA NA NA
## 2009-01-26 1.271524 NA NA NA NA NA
## 2009-02-02 1.289472 NA NA NA NA NAtail(pred.B.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 2.764701 2.965500 2.752886 2.976877 2.636805 2.660384
## 2020-11-09 NA 1.873498 2.017982 1.867302 2.013272 1.776047 1.839694
## 2020-11-16 NA 2.792068 3.047596 2.849108 3.040833 2.705263 2.812665
## 2020-11-23 NA 2.603292 2.838799 2.679337 2.823850 2.522378 2.649897
## 2020-11-30 NA 3.163877 3.473374 3.328649 3.475621 3.126971 3.312110
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2020-11-02 3.062673 3.176020 3.223085 2.786930 3.076671 3.712724 2.774094
## 2020-11-09 2.065239 2.077058 2.150099 1.909664 2.079671 2.342261 1.898749
## 2020-11-16 3.084418 3.100381 3.191465 2.918309 3.099083 3.409805 2.913334
## 2020-11-23 2.830146 2.779059 2.931651 2.750170 2.852951 3.048408 2.734960
## 2020-11-30 3.409974 3.378018 3.541922 3.413220 3.425658 3.574095 3.416408
## 2020-12-07 NA NA NA NA NA NA NA
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2020-11-02 2.877516 3.010112 3.039629 2.719043 2.891563 2.540435 2.811946
## 2020-11-09 1.942484 2.013318 2.039921 1.880071 1.940376 1.728267 1.921446
## 2020-11-16 2.915050 3.001387 3.080285 2.890828 2.963386 2.599209 2.895109
## 2020-11-23 2.708290 2.776693 2.825854 2.732935 2.786921 2.414453 2.707351
## 2020-11-30 3.302229 3.363341 3.445900 3.411289 3.446248 2.955127 3.323752
## 2020-12-07 NA NA NA NA NA NA NA
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2020-11-02 3.053207 2.726051 2.666929 2.625885 2.922180 2.958514 2.661246
## 2020-11-09 2.046592 1.881399 1.802893 1.783246 2.001402 2.033237 1.830403
## 2020-11-16 3.086108 2.886113 2.715707 2.665841 3.046345 3.114578 2.805759
## 2020-11-23 2.849567 2.723671 2.521179 2.486540 2.868559 2.940123 2.649933
## 2020-11-30 3.494329 3.367009 3.086139 3.024817 3.543843 3.664363 3.296532
## 2020-12-07 NA NA NA NA NA NA NA
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2020-11-02 2.634251 2.936375 2.622613 2.836325 2.460262 2.561170 2.773820
## 2020-11-09 1.790993 1.984224 1.836482 1.895818 1.719330 1.774850 1.882119
## 2020-11-16 2.725861 2.992730 2.798335 2.849912 2.548613 2.758144 2.859204
## 2020-11-23 2.544480 2.775065 2.670692 2.685115 2.396096 2.618570 2.711090
## 2020-11-30 3.139794 3.402830 3.317087 3.257333 2.913272 3.301905 3.311188
## 2020-12-07 NA NA NA NA NA NA NA
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2020-11-02 2.692386 2.818715 2.680553 2.624194 2.874322 2.815077 2.768683
## 2020-11-09 1.851665 1.905966 1.860481 1.763521 1.960487 1.937454 1.886314
## 2020-11-16 2.838348 2.879669 2.882865 2.622652 2.937517 2.955529 2.890088
## 2020-11-23 2.673244 2.695895 2.741832 2.423400 2.730532 2.791597 2.700770
## 2020-11-30 3.327800 3.305835 3.450326 2.937237 3.284964 3.462568 3.348806
## 2020-12-07 NA NA NA NA NA NA NA
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 2.854009 2.527001 2.799377 2.776459 3.212378 2.817308 2.696965
## 2020-11-09 1.960736 1.712375 1.939768 1.904529 2.083031 1.923895 1.865104
## 2020-11-16 2.928421 2.552294 2.974432 2.931535 3.110215 2.914900 2.858597
## 2020-11-23 2.750455 2.367249 2.810814 2.775721 2.888619 2.733047 2.735709
## 2020-11-30 3.367352 2.876648 3.476432 3.469285 3.450417 3.370071 3.374917
## 2020-12-07 NA NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 2.841391 2.720054 3.049047 2.712413 2.796957 2.985805 2.881761
## 2020-11-09 1.936326 1.900091 2.053057 1.852071 1.917061 2.020629 1.943521
## 2020-11-16 2.890544 2.913881 3.085088 2.828743 2.901677 3.019545 2.892786
## 2020-11-23 2.710763 2.779251 2.865220 2.662757 2.730022 2.817579 2.694669
## 2020-11-30 3.311958 3.473212 3.489430 3.289792 3.364237 3.443822 3.277358
## 2020-12-07 NA NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2020-11-02 2.858291 2.599410 2.977927 2.746303 2.715763 3.016283 3.007128
## 2020-11-09 1.952332 1.816145 1.983411 1.858735 1.852717 2.012983 2.034449
## 2020-11-16 2.960321 2.812656 3.004052 2.795768 2.820002 3.011793 3.048236
## 2020-11-23 2.778851 2.689736 2.806470 2.598910 2.654363 2.782023 2.842311
## 2020-11-30 3.427676 3.389900 3.462037 3.176563 3.207709 3.376276 3.479447
## 2020-12-07 NA NA NA NA NA NA NA
## d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.843326 NA NA NA NA NA
## 2020-11-09 1.886420 NA NA NA NA NA
## 2020-11-16 2.784749 NA NA NA NA NA
## 2020-11-23 2.543297 NA NA NA NA NA
## 2020-11-30 3.053890 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NAsummary(pred.B.cv)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX
## obs 0.13728827 0.13728827 0.10945441
## average 0.08703924 0.08703924 0.06666056
##
## R2:
## EX.mu EX.mu.beta EX
## obs 0.7201945 0.7201945 0.8221491
## average 0.6634540 0.6634540 0.8025977
##
## Coverage of 95% prediction intervals:
## EX
## obs 0.9417722
## average 0.8730159summary(pred.B.cv.log)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.1926108 0.1926108 0.14780488 0.14784003
## average 0.1032672 0.1032672 0.06923018 0.06983377
##
## R2:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.6970259 0.6970259 0.8215888 0.8215040
## average 0.6830970 0.6830970 0.8575730 0.8550786
##
## Coverage of 95% prediction intervals:
## EX.pred
## obs 0.9417722
## average 0.8730159par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModB.jpeg"),
width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()## quartz_off_screen
## 23.2.5 Predictions for 2009-2020
What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites
x.B <- coef(est.denver.model.B, pars = "cov")$par
x.B## [1] -15.000000 -7.107261 -7.617030 13.498730 -2.979686 -4.934727E.B <- predict(denver.model.B, est.denver.model.B, LTA = T,
transform="unbiased", pred.var=FALSE)
print(E.B)## Prediction for STmodel.
##
## Regression parameters:
## 0 Spatio-temporal covariate(s).
## 39 beta-fields regression parameters in x$pars.
##
## Regression parameters are assumed to be known and
## prediction variances do NOT include
## uncertainties in regression parameters.
##
## Prediction of beta-fields, (x$beta):
## List of 2
## $ mu: num [1:69, 1:3] 0.396 0.199 0.268 0.233 0.247 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX: num [1:69, 1:3] 0.396 0.199 0.268 0.233 0.247 ...
## ..- attr(*, "dimnames")=List of 2
##
## Predictions for log-Gaussian field of type: unbiased
##
## Predictions for 622 times at 69 locations.
## List of 5
## $ EX.mu : num [1:622, 1:69] 1.49 1.52 1.61 1.56 1.57 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.mu.beta: num [1:622, 1:69] 1.53 1.56 1.65 1.59 1.6 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX : num [1:622, 1:69] 1.57 1.6 1.69 1.63 1.64 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.pred : num [1:622, 1:69] 1.58 1.61 1.7 1.64 1.65 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.EX : num [1:622, 1:69] 0.427 0.444 0.501 0.464 0.471 ...
## ..- attr(*, "dimnames")=List of 2
##
## Variances have been computed.
## List of 2
## $ log.VX : num [1:622, 1:69] 0.0513 0.0512 0.0511 0.051 0.051 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.VX.pred: num [1:622, 1:69] 0.0585 0.0584 0.0583 0.0582 0.0582 ...
## ..- attr(*, "dimnames")=List of 2
##
## Mean squared prediciton errors NOT been computed.
##
## 69 temporal averages have been compute.
## List of 1
## $ LTA:'data.frame': 69 obs. of 4 variables:week_preds <- as.data.frame(E.B$EX) %>%
mutate(week = as.Date(rownames(E.B$EX)),
year = year(as.Date(rownames(E.B$EX))))
week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]
box_preds <- week_preds %>%
select(week, all_of(week_sites)) %>%
#filter(week %in% week_weeks) %>%
mutate(month = month(week),
year = year(week)) %>%
filter(year %in% c(2009:2020))
box_data <- box_preds %>%
pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)## week month year location
## Min. :2009-01-05 Min. : 1.000 Min. :2009 Length:42228
## 1st Qu.:2012-01-09 1st Qu.: 4.000 1st Qu.:2012 Class :character
## Median :2014-12-29 Median : 7.000 Median :2014 Mode :character
## Mean :2014-12-28 Mean : 6.494 Mean :2014
## 3rd Qu.:2017-12-18 3rd Qu.: 9.000 3rd Qu.:2017
## Max. :2020-12-07 Max. :12.000 Max. :2020
##
## pred
## Min. :0.3575
## 1st Qu.:1.1777
## Median :1.4767
## Mean :1.5974
## 3rd Qu.:2.0105
## Max. :4.0755
## NA's :861#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
group_by(month) %>%
summarize(median = median(pred, na.rm=T)) %>%
arrange(desc(median))
box_summarypred_box_plot <- ggplot(box_data) +
geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
show.legend = F) +
#scale_color_viridis(name = "Month", discrete = T) +
facet_wrap(. ~ year, scales = "free") +
xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
# theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot
4 Summary Table
| model | trend_data | cov_beta0 | cov_beta1 | cov_beta2 | cov_nu_error | no_basis_fns | df_per_year | all_converge | cv_rmse_obs | cv_rmse_avg | cv_r2_obs | cv_r2_avg | cv_coverage_obs | cv_coverage_avg |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Model A | BC + NO2 + PM2.5 | iid | iid | NA | exp | 1 | 4 | TRUE | 0.15 | 0.06 | 0.81 | 0.89 | 0.94 | 0.89 |
| Model B | BC + NO2 + PM2.5 | iid | iid | iid | exp | 2 | 4 | TRUE | 0.15 | 0.07 | 0.82 | 0.86 | 0.94 | 0.87 |